stanEx2
options(scipen=100,digits=3)
library(cmdstanr)
options(mc.cores=parallel::detectCores())
options(cmdstanr_max_rows=1000)
library(gridExtra)execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(4,1),mar=c(1,5,1,1))
drw[,1,] |> plot(type='l',xlab='',ylab=pr)
drw[,2,] |> plot(type='l',xlab='',ylab=pr)
drw[,3,] |> plot(type='l',xlab='',ylab=pr)
drw[,4,] |> plot(type='l',xlab='',ylab=pr)
par(mar=c(3,5,3,3))
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}regression
normal linear model
ex5-1.stan
normal regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~normal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=normal_rng(m1[i],s);
}
}n=30
mdl=cmdstan_model('./ex5-1.stan')single regression
tb=tibble(x=runif(n,0,9),
y=rnorm(n,x,1))
f=formula(y~x)
par(mfrow=c(1,1))
plot(tb$x,tb$y)
qplot(data=tb,x,y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -115.099
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 20 -13.7106 0.000154046 0.000912282 0.9892 0.9892 24
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.mle## variable estimate
## lp__ -13.71
## b[1] -0.05
## b[2] 0.97
## s 0.96
## y1[1] 8.70
## y1[2] 9.13
## y1[3] 8.41
## y1[4] 6.23
## y1[5] 8.22
## y1[6] 5.17
## y1[7] 3.04
## y1[8] 2.88
## y1[9] 6.78
## y1[10] 6.76
## y1[11] 1.72
## y1[12] 7.52
## y1[13] 6.35
## y1[14] 5.34
## y1[15] 2.95
## y1[16] 3.87
## y1[17] 0.13
## y1[18] 2.68
## y1[19] 1.97
## y1[20] 3.53
## y1[21] 3.88
## y1[22] 4.25
## y1[23] 1.04
## y1[24] -0.71
## y1[25] 0.00
## y1[26] 0.41
## y1[27] 2.23
## y1[28] 5.19
## y1[29] -0.12
## y1[30] 2.67
## m1[1] 8.39
## m1[2] 8.69
## m1[3] 6.75
## m1[4] 7.05
## m1[5] 7.88
## m1[6] 5.92
## m1[7] 3.61
## m1[8] 2.98
## m1[9] 5.59
## m1[10] 8.23
## m1[11] 1.99
## m1[12] 4.84
## m1[13] 5.32
## m1[14] 5.93
## m1[15] 2.87
## m1[16] 3.38
## m1[17] -0.03
## m1[18] 2.60
## m1[19] 3.58
## m1[20] 4.37
## m1[21] 3.77
## m1[22] 5.93
## m1[23] 1.89
## m1[24] 0.17
## m1[25] 0.94
## m1[26] 0.57
## m1[27] 2.28
## m1[28] 4.86
## m1[29] 1.39
## m1[30] 2.46mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.35 -14.99 1.36 1.09 -18.00 -13.90 1.00 697 829
## b[1] -0.05 -0.06 0.39 0.37 -0.69 0.59 1.02 315 494
## b[2] 0.97 0.97 0.08 0.07 0.84 1.10 1.02 392 658
## s 1.04 1.02 0.14 0.13 0.83 1.29 1.00 1013 769
## y1[1] 8.37 8.35 1.13 1.09 6.56 10.20 1.00 1836 1661
## y1[2] 8.69 8.70 1.14 1.10 6.80 10.52 1.00 1979 2013
## y1[3] 6.77 6.76 1.07 1.04 5.04 8.49 1.00 1799 1937
## y1[4] 7.03 7.04 1.08 1.03 5.26 8.81 1.00 2026 1971
## y1[5] 7.88 7.90 1.12 1.07 6.06 9.68 1.00 2066 1807
## y1[6] 5.89 5.91 1.07 1.08 4.08 7.58 1.00 1918 1893
## y1[7] 3.65 3.66 1.08 1.08 1.90 5.38 1.00 2001 1820
## y1[8] 3.00 3.00 1.07 1.03 1.23 4.78 1.00 1686 1721
## y1[9] 5.58 5.59 1.09 1.08 3.75 7.36 1.00 1958 1905
## y1[10] 8.24 8.23 1.16 1.14 6.32 10.16 1.00 2056 1771
## y1[11] 2.01 1.97 1.10 1.07 0.23 3.88 1.00 1643 1711
## y1[12] 4.85 4.86 1.06 1.01 3.14 6.66 1.00 1848 1908
## y1[13] 5.34 5.35 1.07 1.01 3.56 7.08 1.00 1846 1990
## y1[14] 5.92 5.92 1.06 1.00 4.13 7.66 1.00 1881 1855
## y1[15] 2.86 2.86 1.06 1.03 1.17 4.59 1.00 1743 1931
## y1[16] 3.38 3.40 1.10 1.07 1.56 5.19 1.00 2042 1562
## y1[17] -0.05 -0.07 1.15 1.13 -1.85 1.87 1.00 1320 1425
## y1[18] 2.60 2.59 1.07 1.08 0.88 4.41 1.00 1691 1743
## y1[19] 3.62 3.59 1.07 1.01 1.87 5.34 1.00 1813 1757
## y1[20] 4.40 4.40 1.06 1.03 2.62 6.13 1.00 2061 1969
## y1[21] 3.73 3.73 1.05 1.04 2.03 5.48 1.00 2014 1930
## y1[22] 5.96 5.98 1.06 1.04 4.22 7.67 1.00 1988 1998
## y1[23] 1.89 1.86 1.08 1.10 0.08 3.66 1.00 1871 1974
## y1[24] 0.19 0.18 1.09 1.10 -1.51 1.94 1.00 1443 1837
## y1[25] 0.91 0.90 1.11 1.10 -0.85 2.74 1.00 1155 1573
## y1[26] 0.60 0.57 1.09 1.04 -1.25 2.37 1.00 1590 1690
## y1[27] 2.28 2.29 1.06 1.04 0.57 3.98 1.00 1955 1807
## y1[28] 4.84 4.84 1.06 1.00 3.17 6.61 1.00 1913 1853
## y1[29] 1.45 1.45 1.10 1.07 -0.40 3.26 1.00 1650 1646
## y1[30] 2.45 2.47 1.06 1.04 0.69 4.15 1.00 1899 1793
## m1[1] 8.39 8.39 0.39 0.37 7.73 9.05 1.01 1003 1238
## m1[2] 8.69 8.69 0.42 0.38 7.99 9.39 1.01 945 1256
## m1[3] 6.75 6.74 0.29 0.28 6.28 7.23 1.00 1663 1544
## m1[4] 7.05 7.05 0.31 0.29 6.54 7.56 1.00 1469 1479
## m1[5] 7.88 7.88 0.36 0.34 7.29 8.49 1.00 1128 1299
## m1[6] 5.92 5.92 0.25 0.23 5.52 6.33 1.00 2172 1521
## m1[7] 3.61 3.61 0.21 0.20 3.27 3.96 1.01 737 1038
## m1[8] 2.98 2.99 0.22 0.21 2.62 3.35 1.01 501 836
## m1[9] 5.59 5.59 0.23 0.22 5.22 5.98 1.00 2216 1467
## m1[10] 8.23 8.23 0.38 0.36 7.59 8.87 1.00 1038 1290
## m1[11] 1.99 2.00 0.26 0.24 1.55 2.42 1.02 367 647
## m1[12] 4.84 4.84 0.21 0.20 4.51 5.18 1.00 2037 1613
## m1[13] 5.32 5.32 0.22 0.21 4.97 5.69 1.00 2202 1592
## m1[14] 5.93 5.93 0.25 0.23 5.53 6.34 1.00 2171 1521
## m1[15] 2.87 2.88 0.22 0.21 2.50 3.24 1.01 476 781
## m1[16] 3.38 3.38 0.21 0.21 3.03 3.73 1.01 625 934
## m1[17] -0.03 -0.04 0.38 0.37 -0.67 0.61 1.02 315 489
## m1[18] 2.60 2.60 0.23 0.22 2.21 2.99 1.01 430 715
## m1[19] 3.58 3.58 0.21 0.20 3.24 3.92 1.01 718 1027
## m1[20] 4.37 4.37 0.20 0.20 4.05 4.70 1.00 1491 1522
## m1[21] 3.76 3.77 0.20 0.20 3.43 4.10 1.00 832 1110
## m1[22] 5.93 5.93 0.25 0.23 5.52 6.34 1.00 2171 1521
## m1[23] 1.89 1.89 0.27 0.25 1.44 2.33 1.02 361 638
## m1[24] 0.17 0.17 0.37 0.35 -0.45 0.78 1.02 316 479
## m1[25] 0.94 0.94 0.32 0.31 0.40 1.48 1.02 326 549
## m1[26] 0.57 0.57 0.35 0.33 -0.01 1.14 1.02 320 534
## m1[27] 2.28 2.28 0.25 0.23 1.86 2.68 1.02 391 670
## m1[28] 4.86 4.86 0.21 0.21 4.53 5.20 1.00 2046 1612
## m1[29] 1.39 1.39 0.29 0.28 0.89 1.88 1.02 338 579
## m1[30] 2.46 2.46 0.24 0.23 2.06 2.85 1.01 411 715
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
multiple regression
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,x1-x2,1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -36.1196
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -13.6854 0.000193542 0.00129532 1 1 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -13.69
## b[1] 1.40
## b[2] 0.84
## b[3] -1.12
## s 0.96
## y1[1] -1.28
## y1[2] 2.02
## y1[3] -4.52
## y1[4] 5.49
## y1[5] -3.82
## y1[6] -2.10
## y1[7] -1.27
## y1[8] 5.65
## y1[9] -0.65
## y1[10] 2.97
## y1[11] 1.38
## y1[12] 1.75
## y1[13] 8.74
## y1[14] 2.60
## y1[15] -2.98
## y1[16] 3.43
## y1[17] 0.47
## y1[18] 6.46
## y1[19] -0.90
## y1[20] 1.03
## y1[21] 3.67
## y1[22] 1.96
## y1[23] -6.01
## y1[24] 0.90
## y1[25] -0.11
## y1[26] 3.25
## y1[27] -2.32
## y1[28] -1.46
## y1[29] 4.81
## y1[30] 5.89
## m1[1] -1.60
## m1[2] 0.46
## m1[3] -2.91
## m1[4] 5.70
## m1[5] -5.14
## m1[6] -1.86
## m1[7] -1.59
## m1[8] 6.31
## m1[9] -0.22
## m1[10] 3.90
## m1[11] 3.28
## m1[12] 2.74
## m1[13] 7.93
## m1[14] 2.66
## m1[15] -2.73
## m1[16] 2.06
## m1[17] 1.27
## m1[18] 6.81
## m1[19] -0.90
## m1[20] 2.91
## m1[21] 2.36
## m1[22] 2.05
## m1[23] -5.62
## m1[24] -0.12
## m1[25] -1.31
## m1[26] 1.37
## m1[27] -2.45
## m1[28] -1.83
## m1[29] 6.78
## m1[30] 4.90mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -15.90 -15.58 1.49 1.36 -18.91 -14.13 1.01 513 399
## b[1] 1.41 1.44 0.75 0.77 0.22 2.58 1.01 255 262
## b[2] 0.84 0.83 0.09 0.09 0.69 0.99 1.02 306 373
## b[3] -1.12 -1.12 0.09 0.08 -1.27 -0.98 1.00 510 590
## s 1.06 1.05 0.15 0.14 0.85 1.34 1.00 1105 941
## y1[1] -1.60 -1.60 1.13 1.13 -3.51 0.24 1.00 2036 1983
## y1[2] 0.43 0.41 1.10 1.07 -1.37 2.26 1.00 1898 1803
## y1[3] -2.92 -2.93 1.15 1.11 -4.77 -1.01 1.00 1939 2002
## y1[4] 5.71 5.74 1.12 1.08 3.88 7.56 1.00 1933 2014
## y1[5] -5.13 -5.12 1.16 1.13 -7.02 -3.21 1.00 1455 1332
## y1[6] -1.85 -1.86 1.14 1.17 -3.72 -0.04 1.00 1612 1781
## y1[7] -1.63 -1.66 1.10 1.06 -3.44 0.19 1.00 1929 1822
## y1[8] 6.33 6.32 1.17 1.11 4.38 8.21 1.00 1615 1724
## y1[9] -0.20 -0.21 1.15 1.14 -2.07 1.67 1.00 1385 1956
## y1[10] 3.96 3.97 1.10 1.04 2.13 5.81 1.00 1926 1919
## y1[11] 3.28 3.26 1.13 1.08 1.43 5.18 1.00 2023 2097
## y1[12] 2.74 2.72 1.10 1.06 0.96 4.54 1.00 1711 1853
## y1[13] 7.94 7.92 1.16 1.12 6.15 9.86 1.00 2091 1889
## y1[14] 2.70 2.72 1.23 1.16 0.60 4.68 1.00 914 1402
## y1[15] -2.72 -2.71 1.16 1.17 -4.62 -0.79 1.00 1476 1529
## y1[16] 2.11 2.13 1.13 1.11 0.20 3.92 1.00 1711 1702
## y1[17] 1.28 1.30 1.11 1.08 -0.60 3.07 1.00 1739 1915
## y1[18] 6.78 6.81 1.12 1.10 4.96 8.64 1.00 1990 1961
## y1[19] -0.93 -0.93 1.07 1.03 -2.67 0.81 1.00 1955 1738
## y1[20] 2.94 2.91 1.09 1.04 1.12 4.75 1.00 2102 1973
## y1[21] 2.40 2.43 1.14 1.10 0.46 4.25 1.00 1619 1723
## y1[22] 2.05 2.05 1.09 1.05 0.29 3.86 1.00 1978 1689
## y1[23] -5.63 -5.64 1.17 1.14 -7.50 -3.75 1.00 1832 1857
## y1[24] -0.12 -0.11 1.10 1.06 -1.91 1.71 1.00 1880 1930
## y1[25] -1.31 -1.33 1.12 1.07 -3.14 0.53 1.00 1530 1886
## y1[26] 1.37 1.37 1.08 1.06 -0.37 3.16 1.00 2004 1931
## y1[27] -2.48 -2.48 1.11 1.07 -4.36 -0.66 1.00 1737 1932
## y1[28] -1.86 -1.85 1.18 1.15 -3.75 0.10 1.00 912 1227
## y1[29] 6.78 6.75 1.14 1.08 4.89 8.68 1.00 1799 1733
## y1[30] 4.89 4.89 1.14 1.10 3.05 6.76 1.00 1640 1839
## m1[1] -1.61 -1.61 0.27 0.27 -2.05 -1.14 1.00 1472 1389
## m1[2] 0.45 0.45 0.31 0.31 -0.06 0.95 1.00 600 964
## m1[3] -2.92 -2.91 0.34 0.33 -3.49 -2.35 1.00 1152 1412
## m1[4] 5.71 5.70 0.31 0.31 5.20 6.22 1.00 1683 1322
## m1[5] -5.14 -5.13 0.46 0.45 -5.91 -4.38 1.01 506 623
## m1[6] -1.86 -1.87 0.34 0.34 -2.44 -1.30 1.00 748 1297
## m1[7] -1.59 -1.59 0.25 0.24 -2.01 -1.18 1.00 1874 1488
## m1[8] 6.31 6.32 0.41 0.39 5.63 6.98 1.00 615 648
## m1[9] -0.22 -0.23 0.37 0.38 -0.82 0.37 1.01 485 894
## m1[10] 3.90 3.91 0.30 0.29 3.42 4.39 1.00 919 1246
## m1[11] 3.28 3.28 0.23 0.22 2.92 3.65 1.00 1813 1508
## m1[12] 2.74 2.74 0.25 0.24 2.33 3.15 1.01 570 454
## m1[13] 7.93 7.92 0.41 0.40 7.27 8.61 1.00 1684 1267
## m1[14] 2.67 2.69 0.54 0.54 1.80 3.53 1.01 271 279
## m1[15] -2.73 -2.72 0.35 0.34 -3.32 -2.15 1.01 432 422
## m1[16] 2.06 2.04 0.41 0.41 1.39 2.70 1.01 403 701
## m1[17] 1.28 1.28 0.28 0.27 0.83 1.72 1.01 358 344
## m1[18] 6.82 6.80 0.36 0.35 6.25 7.40 1.00 1857 1269
## m1[19] -0.90 -0.91 0.23 0.22 -1.28 -0.53 1.00 1905 1365
## m1[20] 2.91 2.91 0.24 0.23 2.52 3.30 1.00 1706 1260
## m1[21] 2.36 2.37 0.34 0.34 1.80 2.90 1.00 499 985
## m1[22] 2.05 2.06 0.21 0.20 1.71 2.41 1.00 1019 865
## m1[23] -5.62 -5.62 0.42 0.40 -6.31 -4.91 1.00 1216 1388
## m1[24] -0.12 -0.12 0.21 0.20 -0.47 0.23 1.01 1024 814
## m1[25] -1.31 -1.32 0.35 0.35 -1.89 -0.74 1.00 631 1087
## m1[26] 1.37 1.37 0.20 0.20 1.04 1.71 1.00 2062 1444
## m1[27] -2.46 -2.46 0.31 0.30 -2.98 -1.94 1.01 556 603
## m1[28] -1.83 -1.81 0.51 0.51 -2.66 -1.01 1.01 277 277
## m1[29] 6.79 6.79 0.37 0.37 6.17 7.40 1.00 1143 1161
## m1[30] 4.91 4.92 0.38 0.37 4.29 5.54 1.01 443 450
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANOVA
tb=tibble(c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,(c=='b')*2-(c=='c')*2,1))
f=formula(y~c)
par(mfrow=c(1,1))
boxplot(y~c,tb)
qplot(data=tb,c,y,geom='boxplot')
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -52.9384
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 14 -16.6315 0.00024107 0.00123899 0.9856 0.9856 22
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -16.63
## b[1] -0.18
## b[2] 2.15
## b[3] -2.35
## s 1.06
## y1[1] -0.42
## y1[2] 3.05
## y1[3] 0.80
## y1[4] -0.22
## y1[5] -1.48
## y1[6] -0.75
## y1[7] 3.59
## y1[8] 1.46
## y1[9] 2.83
## y1[10] -0.17
## y1[11] 0.86
## y1[12] -1.61
## y1[13] -0.23
## y1[14] -0.40
## y1[15] 1.04
## y1[16] 2.47
## y1[17] 0.66
## y1[18] 0.42
## y1[19] -3.97
## y1[20] -1.22
## y1[21] -1.89
## y1[22] 0.18
## y1[23] -1.63
## y1[24] -3.74
## y1[25] 1.34
## y1[26] 2.07
## y1[27] -0.87
## y1[28] -3.20
## y1[29] 0.99
## y1[30] -3.61
## m1[1] -2.53
## m1[2] 1.97
## m1[3] 1.97
## m1[4] -0.18
## m1[5] -0.18
## m1[6] -2.53
## m1[7] 1.97
## m1[8] 1.97
## m1[9] 1.97
## m1[10] -2.53
## m1[11] -0.18
## m1[12] -2.53
## m1[13] -0.18
## m1[14] -2.53
## m1[15] 1.97
## m1[16] 1.97
## m1[17] 1.97
## m1[18] -0.18
## m1[19] -2.53
## m1[20] -0.18
## m1[21] -0.18
## m1[22] -0.18
## m1[23] -2.53
## m1[24] -2.53
## m1[25] 1.97
## m1[26] 1.97
## m1[27] -2.53
## m1[28] -0.18
## m1[29] -0.18
## m1[30] -2.53mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -18.70 -18.36 1.52 1.34 -21.56 -16.91 1.00 629 943
## b[1] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## b[2] 2.15 2.17 0.52 0.50 1.30 3.02 1.00 712 744
## b[3] -2.33 -2.34 0.52 0.52 -3.19 -1.48 1.01 753 807
## s 1.17 1.15 0.17 0.16 0.93 1.48 1.00 1432 1442
## y1[1] -2.49 -2.45 1.24 1.25 -4.53 -0.50 1.00 1731 1821
## y1[2] 1.98 1.97 1.26 1.21 -0.11 4.02 1.00 1894 2001
## y1[3] 1.94 1.94 1.24 1.20 -0.07 3.95 1.00 1901 1923
## y1[4] -0.15 -0.15 1.25 1.19 -2.20 1.89 1.00 1731 1612
## y1[5] -0.21 -0.19 1.27 1.26 -2.33 1.80 1.00 1691 2101
## y1[6] -2.58 -2.57 1.24 1.21 -4.57 -0.54 1.00 1861 1842
## y1[7] 1.97 1.98 1.27 1.26 -0.09 4.11 1.00 1800 1806
## y1[8] 1.94 1.90 1.22 1.22 -0.02 3.99 1.00 1795 1854
## y1[9] 1.92 1.94 1.21 1.15 -0.03 3.97 1.00 2008 1875
## y1[10] -2.53 -2.53 1.24 1.21 -4.54 -0.53 1.00 1832 1637
## y1[11] -0.20 -0.19 1.22 1.19 -2.21 1.84 1.00 1823 1841
## y1[12] -2.51 -2.52 1.25 1.26 -4.59 -0.48 1.00 2305 1876
## y1[13] -0.20 -0.19 1.23 1.25 -2.22 1.81 1.00 1622 1875
## y1[14] -2.56 -2.55 1.26 1.22 -4.64 -0.49 1.00 2057 1615
## y1[15] 1.99 1.99 1.23 1.20 0.00 4.03 1.00 2121 1841
## y1[16] 1.96 1.96 1.22 1.19 -0.06 3.92 1.00 1903 1990
## y1[17] 1.98 1.99 1.24 1.22 -0.05 4.02 1.00 1828 1984
## y1[18] -0.18 -0.20 1.22 1.23 -2.12 1.81 1.00 1803 1916
## y1[19] -2.53 -2.49 1.26 1.20 -4.59 -0.48 1.00 1830 1915
## y1[20] -0.18 -0.20 1.27 1.23 -2.18 1.90 1.00 1640 1390
## y1[21] -0.18 -0.19 1.26 1.23 -2.19 1.89 1.00 1918 1348
## y1[22] -0.17 -0.20 1.28 1.21 -2.37 1.89 1.00 1754 1819
## y1[23] -2.55 -2.54 1.24 1.21 -4.57 -0.62 1.00 1982 1713
## y1[24] -2.51 -2.51 1.24 1.19 -4.56 -0.51 1.00 1898 1796
## y1[25] 1.95 1.98 1.25 1.22 -0.01 4.02 1.00 1787 1886
## y1[26] 1.99 1.98 1.26 1.18 -0.06 4.04 1.00 2025 1698
## y1[27] -2.52 -2.55 1.25 1.20 -4.50 -0.45 1.00 1854 1934
## y1[28] -0.22 -0.23 1.24 1.18 -2.27 1.82 1.00 1973 1880
## y1[29] -0.20 -0.19 1.26 1.23 -2.28 1.83 1.00 1932 1931
## y1[30] -2.50 -2.48 1.22 1.18 -4.54 -0.54 1.00 1937 1821
## m1[1] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[2] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[3] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[4] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[5] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[6] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[7] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[8] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[9] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[10] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[11] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[12] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[13] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[14] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[15] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[16] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[17] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[18] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[19] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[20] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[21] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[22] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[23] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[24] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[25] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[26] 1.96 1.97 0.37 0.36 1.36 2.57 1.00 1130 1219
## m1[27] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
## m1[28] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[29] -0.19 -0.20 0.37 0.35 -0.78 0.41 1.00 820 1042
## m1[30] -2.52 -2.53 0.37 0.36 -3.14 -1.89 1.01 1121 1346
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,shape=c,size=I(2),xlab='obs.',ylab='prd.')
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
ANCOVA
tb=tibble(x=runif(n,0,9),c=sample(c('a','b','c'),n,replace=T),
y=rnorm(n,x+(c=='b')*2-(c=='c')*2,1))
f=formula(y~x+c)
par(mfrow=c(1,1))
#plot(tb$x1,tb$y,pch=tb$c1)
lv=c(0,1,2)
plot(tb$x,tb$y,pch=lv[factor(tb$c)])
qplot(data=tb,x,y,shape=c,size=I(2))
plot(tb$x,tb$y,col=1+lv[factor(tb$c)])
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)## Initial log joint probability = -451.38
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 43 -9.28787 6.39702e-05 0.000959384 1 1 49
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.29
## b[1] 0.40
## b[2] 0.93
## b[3] 1.62
## b[4] -2.19
## s 0.83
## y1[1] 1.90
## y1[2] 5.40
## y1[3] 2.89
## y1[4] -1.32
## y1[5] 0.37
## y1[6] 7.28
## y1[7] -1.70
## y1[8] 9.26
## y1[9] 0.85
## y1[10] 7.72
## y1[11] 4.09
## y1[12] 8.41
## y1[13] 8.72
## y1[14] 3.41
## y1[15] 7.33
## y1[16] 6.15
## y1[17] 6.99
## y1[18] 1.88
## y1[19] 3.16
## y1[20] 3.86
## y1[21] 4.37
## y1[22] 2.42
## y1[23] 7.04
## y1[24] 3.61
## y1[25] 4.29
## y1[26] -0.75
## y1[27] 0.90
## y1[28] 6.71
## y1[29] 6.05
## y1[30] 8.05
## m1[1] 3.51
## m1[2] 4.76
## m1[3] 3.82
## m1[4] -0.35
## m1[5] 0.86
## m1[6] 7.81
## m1[7] -1.12
## m1[8] 7.89
## m1[9] -0.10
## m1[10] 6.65
## m1[11] 3.19
## m1[12] 8.32
## m1[13] 9.06
## m1[14] 5.43
## m1[15] 6.77
## m1[16] 5.75
## m1[17] 7.63
## m1[18] 2.18
## m1[19] 2.84
## m1[20] 4.01
## m1[21] 5.32
## m1[22] 2.75
## m1[23] 6.85
## m1[24] 4.09
## m1[25] 5.75
## m1[26] 0.37
## m1[27] 1.99
## m1[28] 5.94
## m1[29] 6.90
## m1[30] 6.51mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.08 -11.75 1.65 1.44 -15.28 -10.05 1.00 622 840
## b[1] 0.38 0.37 0.44 0.43 -0.31 1.12 1.01 607 702
## b[2] 0.93 0.93 0.07 0.07 0.82 1.04 1.00 1148 1151
## b[3] 1.63 1.62 0.40 0.38 0.98 2.29 1.01 683 961
## b[4] -2.20 -2.19 0.40 0.39 -2.88 -1.55 1.00 643 767
## s 0.93 0.92 0.14 0.13 0.74 1.16 1.00 1868 1480
## y1[1] 3.49 3.47 1.00 1.01 1.87 5.14 1.00 1762 2074
## y1[2] 4.76 4.74 0.97 0.92 3.24 6.38 1.00 1845 2014
## y1[3] 3.80 3.81 1.01 0.99 2.19 5.46 1.00 2019 1922
## y1[4] -0.35 -0.37 1.01 0.98 -2.00 1.30 1.00 1975 1971
## y1[5] 0.89 0.90 1.04 1.04 -0.85 2.56 1.01 1270 1666
## y1[6] 7.80 7.79 0.97 0.98 6.23 9.41 1.00 1878 1846
## y1[7] -1.14 -1.14 1.04 1.04 -2.88 0.52 1.00 1945 1788
## y1[8] 7.89 7.87 1.01 0.98 6.22 9.57 1.00 1792 1917
## y1[9] -0.16 -0.16 0.99 0.95 -1.81 1.46 1.00 1830 1911
## y1[10] 6.67 6.68 1.00 0.96 5.03 8.31 1.00 2054 1959
## y1[11] 3.18 3.19 0.99 0.97 1.55 4.84 1.00 2005 1857
## y1[12] 8.33 8.32 1.01 1.01 6.72 10.00 1.00 2094 1841
## y1[13] 9.09 9.12 1.00 0.95 7.38 10.72 1.00 1775 1559
## y1[14] 5.44 5.44 0.97 1.00 3.94 7.03 1.00 1713 2052
## y1[15] 6.79 6.80 1.01 0.97 5.16 8.42 1.00 2073 1972
## y1[16] 5.76 5.81 0.98 0.96 4.15 7.33 1.00 1994 1930
## y1[17] 7.66 7.64 1.00 0.95 6.07 9.33 1.00 1869 1736
## y1[18] 2.18 2.16 1.04 1.01 0.54 3.91 1.00 1843 1630
## y1[19] 2.85 2.83 0.99 0.93 1.27 4.47 1.00 2051 1802
## y1[20] 4.01 3.98 1.01 0.98 2.44 5.67 1.00 1743 1618
## y1[21] 5.31 5.31 0.99 0.95 3.67 6.91 1.00 1851 1858
## y1[22] 2.71 2.71 0.98 0.99 1.12 4.32 1.00 1679 1879
## y1[23] 6.86 6.83 0.99 0.96 5.18 8.51 1.00 2045 1800
## y1[24] 4.05 4.03 1.00 0.98 2.46 5.68 1.00 1981 2014
## y1[25] 5.75 5.77 0.99 1.00 4.12 7.33 1.00 1930 1721
## y1[26] 0.37 0.35 1.00 0.98 -1.23 2.04 1.00 1969 1866
## y1[27] 1.99 1.98 0.98 0.99 0.43 3.58 1.00 1910 1929
## y1[28] 5.93 5.95 0.96 0.92 4.31 7.45 1.00 1921 1751
## y1[29] 6.90 6.89 0.98 0.97 5.29 8.48 1.00 1988 1892
## y1[30] 6.50 6.50 1.02 1.04 4.78 8.15 1.00 1544 1706
## m1[1] 3.50 3.50 0.30 0.29 3.00 3.99 1.00 1251 1444
## m1[2] 4.76 4.75 0.29 0.28 4.31 5.25 1.00 684 761
## m1[3] 3.81 3.82 0.31 0.30 3.31 4.31 1.00 1232 1489
## m1[4] -0.37 -0.38 0.36 0.35 -0.95 0.20 1.00 1324 1068
## m1[5] 0.84 0.83 0.42 0.40 0.19 1.55 1.01 595 662
## m1[6] 7.81 7.80 0.35 0.34 7.25 8.39 1.00 1415 1426
## m1[7] -1.14 -1.15 0.39 0.38 -1.77 -0.51 1.00 1308 1014
## m1[8] 7.90 7.91 0.33 0.32 7.36 8.42 1.01 1266 1436
## m1[9] -0.12 -0.13 0.35 0.33 -0.68 0.44 1.00 1333 1149
## m1[10] 6.66 6.66 0.30 0.29 6.17 7.15 1.00 1367 1485
## m1[11] 3.19 3.19 0.35 0.33 2.64 3.78 1.00 1383 1294
## m1[12] 8.32 8.32 0.37 0.36 7.73 8.93 1.00 1550 1458
## m1[13] 9.07 9.08 0.37 0.36 8.45 9.67 1.01 1200 1429
## m1[14] 5.42 5.42 0.29 0.28 4.97 5.91 1.00 781 894
## m1[15] 6.77 6.78 0.30 0.29 6.28 7.27 1.00 1358 1386
## m1[16] 5.75 5.74 0.29 0.28 5.27 6.23 1.00 1402 1453
## m1[17] 7.64 7.65 0.32 0.31 7.11 8.16 1.01 1285 1448
## m1[18] 2.17 2.17 0.39 0.38 1.54 2.84 1.00 1354 1196
## m1[19] 2.83 2.83 0.36 0.35 2.26 3.44 1.00 1376 1181
## m1[20] 4.00 3.99 0.30 0.29 3.54 4.51 1.01 615 706
## m1[21] 5.31 5.32 0.36 0.35 4.71 5.90 1.00 1193 1401
## m1[22] 2.74 2.73 0.33 0.32 2.21 3.30 1.01 575 674
## m1[23] 6.85 6.84 0.31 0.31 6.35 7.37 1.00 1143 1286
## m1[24] 4.08 4.09 0.32 0.31 3.57 4.59 1.00 1223 1469
## m1[25] 5.75 5.75 0.29 0.28 5.28 6.24 1.00 1402 1453
## m1[26] 0.34 0.33 0.33 0.32 -0.20 0.89 1.00 1344 1164
## m1[27] 1.97 1.98 0.30 0.29 1.49 2.45 1.00 1318 1243
## m1[28] 5.93 5.93 0.29 0.29 5.47 6.43 1.00 884 946
## m1[29] 6.90 6.89 0.31 0.31 6.40 7.42 1.00 1157 1360
## m1[30] 6.50 6.52 0.41 0.40 5.81 7.18 1.00 1163 1264
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
lv=c(0,1,2)
par(mfrow=c(1,1))
plot(tb$y,tb$y1,pch=lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
plot(tb$y,tb$y1,col=1+lv[factor(tb$c)],xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,col=c,xlab='obs.',ylab='prd.')+
geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black',pch=lv[factor(tb$c)])
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(data=tb,1:n,y,col=c)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
interaction of variable
n=50
mdl=cmdstan_model('./ex5-1.stan')
tb=tibble(x=runif(n,-3,3),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,x+(ca=='b')+(cb=='b')-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,x,y,col=ca),
qplot(data=tb,x,y,col=cb),ncol=2)
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~x+ca)
f1=formula(y~x+ca+cb)
f2=formula(y~x+ca*cb)
fn(f0)## Initial log joint probability = -119.88
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -27.5817 0.000129179 0.000855749 1 1 26
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.64 -29.28 1.52 1.28 -32.48 -27.88 1.01 908 1002
## b[1] 0.71 0.71 0.22 0.21 0.34 1.09 1.00 586 987
## b[2] 0.97 0.97 0.09 0.09 0.81 1.12 1.00 2127 1383
## b[3] 0.51 0.52 0.32 0.33 -0.01 1.03 1.00 515 861
## s 1.11 1.11 0.12 0.12 0.93 1.33 1.00 1529 1448
## y1[1] -1.71 -1.75 1.17 1.15 -3.59 0.27 1.00 1932 2031
## y1[2] -1.79 -1.78 1.16 1.08 -3.71 0.13 1.00 1943 1966
## y1[3] 0.49 0.51 1.16 1.18 -1.44 2.34 1.00 1921 1933
## y1[4] 1.18 1.16 1.12 1.09 -0.64 3.06 1.00 2034 1929
## y1[5] 1.19 1.20 1.14 1.10 -0.72 3.02 1.00 2003 1896
## y1[6] -0.85 -0.86 1.16 1.11 -2.82 1.09 1.00 1912 1954
## y1[7] 1.02 1.01 1.12 1.11 -0.79 2.87 1.00 1765 1909
## y1[8] 1.24 1.24 1.14 1.09 -0.62 3.14 1.00 2074 1958
## y1[9] 3.20 3.15 1.15 1.17 1.34 5.14 1.00 1808 1850
## y1[10] 2.52 2.54 1.14 1.12 0.60 4.40 1.00 1901 1928
## y1[11] 3.44 3.47 1.17 1.12 1.42 5.32 1.00 2154 1972
## y1[12] 2.97 3.00 1.19 1.21 0.99 4.87 1.00 1419 1862
## y1[13] -1.15 -1.13 1.15 1.15 -3.12 0.67 1.00 2113 2096
## y1[14] 2.75 2.73 1.17 1.12 0.84 4.66 1.00 1741 1785
## y1[15] 3.46 3.47 1.19 1.18 1.50 5.39 1.00 2056 1894
## y1[16] -1.77 -1.78 1.13 1.10 -3.60 0.15 1.00 1818 1734
## y1[17] -1.00 -0.99 1.13 1.07 -2.86 0.79 1.00 1922 1823
## y1[18] 1.58 1.57 1.16 1.17 -0.30 3.45 1.00 2202 2020
## y1[19] 0.37 0.36 1.17 1.13 -1.57 2.34 1.00 1797 1908
## y1[20] 0.59 0.59 1.17 1.19 -1.36 2.45 1.00 1830 1763
## y1[21] 2.92 2.88 1.15 1.16 1.09 4.85 1.00 2047 1966
## y1[22] -1.38 -1.37 1.15 1.15 -3.21 0.51 1.00 1734 1811
## y1[23] 0.78 0.79 1.13 1.10 -1.10 2.60 1.00 1638 1871
## y1[24] -2.03 -2.02 1.16 1.18 -3.93 -0.18 1.00 1958 1859
## y1[25] -1.47 -1.44 1.14 1.16 -3.38 0.41 1.00 2068 1934
## y1[26] 3.38 3.40 1.18 1.19 1.38 5.33 1.00 1901 1835
## y1[27] 0.40 0.40 1.14 1.10 -1.47 2.27 1.00 1832 1932
## y1[28] 1.78 1.78 1.12 1.11 -0.06 3.60 1.00 1997 2024
## y1[29] -0.67 -0.67 1.16 1.14 -2.55 1.23 1.00 1801 1882
## y1[30] -1.75 -1.75 1.16 1.14 -3.69 0.14 1.00 1924 2004
## y1[31] 0.99 1.00 1.13 1.10 -0.91 2.90 1.00 2056 1991
## y1[32] 0.10 0.12 1.15 1.11 -1.83 1.95 1.00 2013 1839
## y1[33] 0.59 0.61 1.15 1.16 -1.26 2.46 1.00 1846 1807
## y1[34] -1.33 -1.37 1.19 1.18 -3.26 0.68 1.00 1918 1929
## y1[35] 3.23 3.21 1.16 1.13 1.37 5.18 1.00 1996 1886
## y1[36] 2.97 2.96 1.15 1.12 1.16 4.93 1.00 1925 1762
## y1[37] -1.10 -1.12 1.20 1.21 -3.05 0.95 1.00 1940 1968
## y1[38] 0.14 0.13 1.16 1.14 -1.71 2.06 1.00 1879 1947
## y1[39] 0.56 0.59 1.16 1.16 -1.43 2.48 1.00 2077 1884
## y1[40] 0.16 0.14 1.12 1.11 -1.63 2.00 1.00 1842 1972
## y1[41] 1.19 1.20 1.14 1.10 -0.72 3.01 1.00 1972 1878
## y1[42] 2.63 2.63 1.14 1.14 0.75 4.54 1.00 1985 2012
## y1[43] -1.29 -1.29 1.13 1.13 -3.10 0.55 1.00 2029 2053
## y1[44] 0.92 0.93 1.14 1.12 -0.98 2.70 1.00 1961 1970
## y1[45] -1.20 -1.20 1.15 1.12 -3.07 0.69 1.00 1749 1754
## y1[46] 1.16 1.19 1.15 1.13 -0.74 2.97 1.00 1931 1921
## y1[47] -2.15 -2.18 1.17 1.15 -4.03 -0.26 1.00 1829 2034
## y1[48] 3.41 3.41 1.17 1.13 1.50 5.35 1.00 1860 1959
## y1[49] 3.49 3.51 1.22 1.19 1.38 5.44 1.00 1765 1886
## y1[50] 1.17 1.19 1.15 1.14 -0.74 3.05 1.00 1920 2010
## m1[1] -1.69 -1.69 0.28 0.27 -2.15 -1.23 1.00 1455 1756
## m1[2] -1.83 -1.83 0.29 0.28 -2.30 -1.35 1.00 1523 1757
## m1[3] 0.48 0.48 0.24 0.23 0.08 0.87 1.00 987 1221
## m1[4] 1.15 1.15 0.23 0.22 0.76 1.52 1.00 1050 1091
## m1[5] 1.16 1.16 0.23 0.23 0.77 1.55 1.00 564 947
## m1[6] -0.83 -0.83 0.24 0.24 -1.23 -0.44 1.00 991 1298
## m1[7] 1.05 1.05 0.23 0.22 0.67 1.42 1.00 1035 1104
## m1[8] 1.24 1.24 0.23 0.22 0.86 1.62 1.00 1060 1110
## m1[9] 3.19 3.19 0.30 0.28 2.69 3.68 1.00 1704 1480
## m1[10] 2.58 2.58 0.26 0.25 2.13 3.01 1.00 1466 1308
## m1[11] 3.46 3.47 0.31 0.30 2.93 3.98 1.00 1808 1503
## m1[12] 3.00 3.01 0.35 0.34 2.43 3.56 1.00 706 1023
## m1[13] -1.13 -1.14 0.25 0.25 -1.55 -0.72 1.00 1172 1592
## m1[14] 2.73 2.74 0.33 0.32 2.19 3.27 1.00 677 1043
## m1[15] 3.40 3.40 0.31 0.30 2.87 3.91 1.00 1784 1529
## m1[16] -1.81 -1.81 0.29 0.28 -2.28 -1.33 1.00 1517 1757
## m1[17] -1.01 -1.01 0.25 0.24 -1.41 -0.61 1.00 1091 1539
## m1[18] 1.59 1.59 0.23 0.22 1.21 1.97 1.00 1129 1134
## m1[19] 0.33 0.33 0.24 0.24 -0.07 0.72 1.00 983 1322
## m1[20] 0.60 0.60 0.22 0.21 0.23 0.97 1.00 601 1008
## m1[21] 2.96 2.96 0.28 0.27 2.48 3.42 1.00 1610 1414
## m1[22] -1.36 -1.36 0.34 0.33 -1.91 -0.82 1.00 1116 1530
## m1[23] 0.74 0.74 0.22 0.21 0.37 1.12 1.00 583 987
## m1[24] -2.00 -2.00 0.30 0.29 -2.49 -1.50 1.00 1607 1663
## m1[25] -1.47 -1.47 0.27 0.26 -1.91 -1.03 1.00 1346 1621
## m1[26] 3.35 3.36 0.37 0.37 2.73 3.96 1.00 739 911
## m1[27] 0.35 0.35 0.24 0.24 -0.06 0.73 1.00 984 1322
## m1[28] 1.79 1.79 0.23 0.23 1.41 2.18 1.00 1186 1175
## m1[29] -0.65 -0.65 0.29 0.28 -1.14 -0.20 1.00 1024 1555
## m1[30] -1.75 -1.75 0.29 0.27 -2.22 -1.28 1.00 1487 1676
## m1[31] 1.03 1.04 0.23 0.22 0.65 1.41 1.00 1033 1103
## m1[32] 0.14 0.14 0.25 0.24 -0.28 0.54 1.00 981 1350
## m1[33] 0.62 0.62 0.22 0.22 0.25 0.99 1.00 598 1008
## m1[34] -1.31 -1.31 0.33 0.32 -1.85 -0.77 1.00 1108 1482
## m1[35] 3.22 3.22 0.30 0.28 2.71 3.71 1.00 1717 1414
## m1[36] 2.96 2.96 0.28 0.27 2.48 3.42 1.00 1609 1414
## m1[37] -1.14 -1.14 0.32 0.31 -1.66 -0.62 1.00 1086 1512
## m1[38] 0.13 0.13 0.22 0.22 -0.23 0.49 1.00 681 1017
## m1[39] 0.55 0.55 0.24 0.23 0.16 0.93 1.00 992 1221
## m1[40] 0.13 0.13 0.22 0.22 -0.23 0.49 1.00 681 1017
## m1[41] 1.23 1.23 0.24 0.24 0.84 1.62 1.00 564 834
## m1[42] 2.64 2.64 0.27 0.26 2.19 3.08 1.00 1490 1399
## m1[43] -1.29 -1.29 0.26 0.25 -1.72 -0.86 1.00 1259 1570
## m1[44] 0.91 0.91 0.23 0.22 0.53 1.29 1.00 572 989
## m1[45] -1.20 -1.20 0.25 0.25 -1.62 -0.78 1.00 1217 1445
## m1[46] 1.19 1.19 0.23 0.22 0.81 1.57 1.00 1054 1117
## m1[47] -2.13 -2.12 0.31 0.30 -2.64 -1.62 1.00 1671 1634
## m1[48] 3.43 3.43 0.38 0.37 2.80 4.04 1.00 748 907
## m1[49] 3.45 3.45 0.38 0.37 2.82 4.06 1.00 750 907
## m1[50] 1.18 1.18 0.23 0.22 0.80 1.56 1.00 1053 1117
fn(f1)## Initial log joint probability = -1090.15
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -26.7666 0.000133064 0.0011364 1 1 20
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -29.37 -29.02 1.67 1.42 -32.63 -27.37 1.00 866 1174
## b[1] 0.47 0.46 0.30 0.30 -0.03 0.95 1.01 529 738
## b[2] 0.98 0.98 0.09 0.09 0.83 1.13 1.00 2201 1535
## b[3] 0.50 0.52 0.33 0.32 -0.06 1.04 1.00 675 844
## b[4] 0.40 0.41 0.31 0.32 -0.11 0.90 1.01 695 816
## s 1.11 1.10 0.12 0.12 0.93 1.33 1.00 2074 1631
## y1[1] -1.54 -1.53 1.17 1.13 -3.44 0.35 1.00 1817 1911
## y1[2] -1.67 -1.71 1.13 1.09 -3.43 0.21 1.00 2168 1897
## y1[3] 0.23 0.26 1.16 1.14 -1.65 2.13 1.00 1931 1933
## y1[4] 1.30 1.29 1.14 1.15 -0.55 3.18 1.00 2087 1927
## y1[5] 1.33 1.34 1.16 1.14 -0.58 3.20 1.00 1757 1669
## y1[6] -0.69 -0.70 1.17 1.17 -2.62 1.25 1.00 1925 1910
## y1[7] 1.21 1.19 1.12 1.11 -0.60 2.99 1.00 1734 1923
## y1[8] 0.97 0.99 1.17 1.15 -0.99 2.93 1.00 1611 1928
## y1[9] 3.31 3.28 1.17 1.12 1.38 5.21 1.00 1884 1791
## y1[10] 2.74 2.74 1.17 1.11 0.86 4.69 1.00 1847 1825
## y1[11] 3.62 3.61 1.16 1.13 1.67 5.57 1.00 1980 1917
## y1[12] 3.24 3.27 1.16 1.15 1.32 5.16 1.00 1857 1833
## y1[13] -0.97 -1.01 1.17 1.13 -2.89 0.95 1.00 2076 2150
## y1[14] 2.49 2.51 1.17 1.13 0.54 4.38 1.00 1913 2011
## y1[15] 3.19 3.22 1.21 1.18 1.15 5.12 1.00 1884 1971
## y1[16] -1.66 -1.67 1.18 1.22 -3.61 0.28 1.00 1997 1927
## y1[17] -0.87 -0.85 1.16 1.16 -2.80 1.03 1.00 2010 2012
## y1[18] 1.37 1.38 1.16 1.19 -0.52 3.31 1.00 2067 1973
## y1[19] 0.01 0.05 1.15 1.14 -1.88 1.85 1.00 1700 1769
## y1[20] 0.79 0.77 1.13 1.11 -1.06 2.68 1.00 1715 1933
## y1[21] 2.69 2.69 1.16 1.18 0.79 4.61 1.00 2086 1821
## y1[22] -1.62 -1.62 1.20 1.18 -3.60 0.40 1.00 1995 1803
## y1[23] 0.48 0.49 1.16 1.13 -1.47 2.38 1.01 1547 1858
## y1[24] -2.22 -2.23 1.19 1.15 -4.17 -0.29 1.00 1727 1747
## y1[25] -1.32 -1.31 1.19 1.17 -3.34 0.57 1.00 1910 2016
## y1[26] 3.11 3.11 1.19 1.18 1.19 5.09 1.00 1728 1617
## y1[27] 0.51 0.50 1.16 1.11 -1.41 2.36 1.00 1914 1932
## y1[28] 1.98 1.97 1.15 1.14 0.09 3.88 1.00 1921 1855
## y1[29] -0.49 -0.50 1.19 1.19 -2.43 1.52 1.00 1859 1639
## y1[30] -2.05 -2.03 1.17 1.15 -3.95 -0.11 1.00 1969 2045
## y1[31] 1.16 1.15 1.16 1.14 -0.72 3.04 1.00 1499 1636
## y1[32] 0.28 0.27 1.14 1.10 -1.63 2.13 1.00 1706 1788
## y1[33] 0.74 0.73 1.13 1.12 -1.04 2.63 1.00 1732 1805
## y1[34] -1.22 -1.19 1.14 1.09 -3.10 0.64 1.00 2087 1924
## y1[35] 3.37 3.38 1.15 1.12 1.51 5.23 1.00 1947 1853
## y1[36] 3.09 3.09 1.12 1.12 1.23 4.86 1.00 2113 2102
## y1[37] -1.04 -1.03 1.17 1.14 -3.04 0.85 1.00 1773 1886
## y1[38] 0.25 0.24 1.15 1.11 -1.59 2.09 1.00 1941 1890
## y1[39] 0.65 0.67 1.13 1.09 -1.19 2.51 1.00 1929 1820
## y1[40] 0.28 0.26 1.14 1.17 -1.61 2.13 1.00 2086 2052
## y1[41] 1.37 1.36 1.12 1.07 -0.45 3.31 1.00 1649 1801
## y1[42] 2.42 2.42 1.14 1.17 0.54 4.27 1.00 1997 1654
## y1[43] -1.54 -1.51 1.16 1.14 -3.51 0.30 1.00 1586 1827
## y1[44] 0.66 0.65 1.13 1.14 -1.18 2.52 1.00 1665 1723
## y1[45] -1.08 -1.06 1.15 1.12 -3.02 0.81 1.00 1975 2035
## y1[46] 1.34 1.34 1.15 1.13 -0.57 3.25 1.00 1901 2058
## y1[47] -2.41 -2.39 1.19 1.11 -4.41 -0.42 1.00 1815 1752
## y1[48] 3.19 3.16 1.16 1.12 1.32 5.09 1.00 1837 1836
## y1[49] 3.66 3.64 1.20 1.20 1.74 5.66 1.00 1911 1892
## y1[50] 0.90 0.88 1.15 1.15 -0.97 2.83 1.00 1804 2016
## m1[1] -1.56 -1.56 0.31 0.30 -2.05 -1.05 1.00 1968 1455
## m1[2] -1.70 -1.70 0.32 0.30 -2.21 -1.18 1.00 2025 1428
## m1[3] 0.22 0.21 0.31 0.30 -0.30 0.72 1.00 959 1209
## m1[4] 1.29 1.30 0.26 0.26 0.86 1.71 1.00 1266 1305
## m1[5] 1.32 1.33 0.28 0.28 0.86 1.77 1.00 1022 1184
## m1[6] -0.69 -0.69 0.27 0.26 -1.13 -0.24 1.00 1613 1522
## m1[7] 1.19 1.20 0.26 0.26 0.77 1.61 1.00 1253 1361
## m1[8] 0.98 0.98 0.30 0.29 0.49 1.47 1.00 956 1156
## m1[9] 3.36 3.36 0.32 0.30 2.83 3.90 1.00 1757 1734
## m1[10] 2.74 2.74 0.29 0.28 2.25 3.22 1.00 1573 1484
## m1[11] 3.63 3.63 0.33 0.33 3.09 4.21 1.00 1829 1587
## m1[12] 3.19 3.18 0.38 0.38 2.57 3.81 1.00 1217 1554
## m1[13] -1.00 -1.00 0.28 0.28 -1.45 -0.53 1.00 1739 1524
## m1[14] 2.51 2.51 0.38 0.38 1.90 3.14 1.00 687 881
## m1[15] 3.16 3.16 0.35 0.34 2.58 3.74 1.00 1291 1468
## m1[16] -1.68 -1.68 0.32 0.30 -2.19 -1.16 1.00 2018 1407
## m1[17] -0.87 -0.87 0.28 0.27 -1.32 -0.41 1.00 1684 1610
## m1[18] 1.34 1.34 0.30 0.29 0.86 1.83 1.00 973 1124
## m1[19] 0.07 0.06 0.31 0.31 -0.45 0.58 1.00 964 1176
## m1[20] 0.75 0.76 0.26 0.26 0.31 1.17 1.00 1085 1239
## m1[21] 2.72 2.72 0.33 0.32 2.18 3.26 1.00 1181 1488
## m1[22] -1.65 -1.65 0.40 0.39 -2.31 -1.00 1.01 1144 1541
## m1[23] 0.50 0.49 0.30 0.30 0.00 0.99 1.01 529 738
## m1[24] -2.27 -2.28 0.37 0.38 -2.86 -1.69 1.01 863 1073
## m1[25] -1.34 -1.34 0.30 0.29 -1.81 -0.83 1.00 1878 1436
## m1[26] 3.14 3.13 0.42 0.41 2.45 3.84 1.00 773 889
## m1[27] 0.48 0.49 0.27 0.27 0.03 0.92 1.00 1225 1256
## m1[28] 1.95 1.94 0.26 0.26 1.50 2.38 1.00 1368 1366
## m1[29] -0.53 -0.52 0.32 0.31 -1.05 -0.02 1.00 1286 1306
## m1[30] -2.03 -2.03 0.36 0.36 -2.59 -1.45 1.01 808 1088
## m1[31] 1.18 1.18 0.26 0.26 0.75 1.60 1.00 1252 1307
## m1[32] 0.27 0.28 0.28 0.28 -0.19 0.72 1.00 1230 1270
## m1[33] 0.77 0.78 0.26 0.27 0.33 1.19 1.00 1081 1256
## m1[34] -1.19 -1.19 0.36 0.35 -1.79 -0.61 1.00 1362 1290
## m1[35] 3.39 3.39 0.32 0.31 2.85 3.93 1.00 1766 1734
## m1[36] 3.13 3.12 0.30 0.29 2.61 3.64 1.00 1687 1700
## m1[37] -1.02 -1.01 0.35 0.34 -1.59 -0.46 1.00 1338 1272
## m1[38] 0.28 0.29 0.26 0.26 -0.14 0.70 1.00 1197 1365
## m1[39] 0.69 0.70 0.27 0.26 0.24 1.12 1.00 1222 1282
## m1[40] 0.28 0.29 0.26 0.26 -0.14 0.70 1.00 1197 1365
## m1[41] 1.39 1.40 0.28 0.28 0.92 1.84 1.00 1021 1194
## m1[42] 2.40 2.40 0.32 0.30 1.89 2.92 1.00 1115 1541
## m1[43] -1.56 -1.56 0.33 0.33 -2.10 -1.01 1.01 715 951
## m1[44] 0.67 0.66 0.31 0.30 0.17 1.16 1.01 532 709
## m1[45] -1.06 -1.06 0.29 0.28 -1.52 -0.58 1.00 1766 1524
## m1[46] 1.33 1.34 0.26 0.25 0.90 1.75 1.00 1272 1305
## m1[47] -2.40 -2.40 0.38 0.38 -3.01 -1.79 1.01 894 1088
## m1[48] 3.22 3.21 0.43 0.42 2.52 3.92 1.00 785 889
## m1[49] 3.64 3.64 0.41 0.41 2.98 4.32 1.00 1283 1610
## m1[50] 0.93 0.93 0.30 0.29 0.44 1.41 1.00 955 1104
fn(f2)## Initial log joint probability = -81.3509
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 29 -25.1811 0.000169817 0.0004775 0.9558 0.9558 36
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -28.35 -27.98 1.87 1.60 -32.07 -26.03 1.00 697 849
## b[1] 0.16 0.15 0.37 0.36 -0.43 0.74 1.00 539 828
## b[2] 0.99 0.99 0.09 0.09 0.84 1.14 1.00 2080 1643
## b[3] 1.16 1.17 0.52 0.52 0.26 1.99 1.00 454 691
## b[4] 0.90 0.91 0.45 0.45 0.15 1.65 1.00 479 899
## b[5] -1.03 -1.04 0.65 0.66 -2.10 0.06 1.01 412 659
## s 1.09 1.08 0.12 0.11 0.92 1.30 1.00 1812 1442
## y1[1] -1.37 -1.36 1.15 1.16 -3.24 0.50 1.00 2018 1968
## y1[2] -1.52 -1.51 1.16 1.14 -3.44 0.35 1.00 1948 1957
## y1[3] 0.54 0.54 1.14 1.14 -1.39 2.42 1.00 1859 1932
## y1[4] 1.08 1.10 1.15 1.11 -0.87 2.96 1.00 2028 2016
## y1[5] 1.47 1.44 1.14 1.13 -0.40 3.32 1.00 1930 1870
## y1[6] -0.56 -0.57 1.09 1.08 -2.35 1.21 1.00 2019 1971
## y1[7] 0.98 1.00 1.15 1.12 -0.94 2.79 1.00 2025 1915
## y1[8] 1.32 1.33 1.16 1.13 -0.58 3.29 1.00 1735 1822
## y1[9] 3.19 3.22 1.14 1.14 1.31 5.02 1.00 2010 1996
## y1[10] 2.65 2.66 1.12 1.08 0.85 4.48 1.00 2023 1932
## y1[11] 3.50 3.52 1.14 1.13 1.68 5.35 1.00 1860 1596
## y1[12] 3.39 3.39 1.14 1.11 1.52 5.27 1.00 1764 1647
## y1[13] -0.86 -0.85 1.12 1.08 -2.72 1.03 1.00 2023 1857
## y1[14] 2.19 2.20 1.15 1.12 0.32 4.12 1.00 1308 1768
## y1[15] 3.53 3.55 1.17 1.17 1.64 5.44 1.00 1779 1850
## y1[16] -1.53 -1.52 1.16 1.14 -3.39 0.39 1.00 1927 1857
## y1[17] -0.73 -0.70 1.12 1.11 -2.63 1.10 1.00 1965 1848
## y1[18] 1.68 1.64 1.17 1.11 -0.29 3.61 1.00 1764 1837
## y1[19] 0.38 0.40 1.16 1.12 -1.59 2.21 1.00 1798 1823
## y1[20] 0.89 0.92 1.13 1.12 -0.95 2.74 1.00 1932 1916
## y1[21] 3.08 3.08 1.17 1.09 1.18 5.05 1.00 1890 1660
## y1[22] -1.36 -1.36 1.16 1.16 -3.20 0.59 1.00 1911 1916
## y1[23] 0.20 0.22 1.14 1.15 -1.69 2.07 1.00 1596 1711
## y1[24] -2.56 -2.58 1.16 1.14 -4.52 -0.67 1.00 1621 1843
## y1[25] -1.20 -1.20 1.13 1.15 -3.04 0.66 1.00 2010 1780
## y1[26] 2.83 2.86 1.22 1.14 0.77 4.79 1.00 1932 1834
## y1[27] 0.28 0.27 1.13 1.11 -1.55 2.13 1.00 2092 1960
## y1[28] 1.78 1.76 1.12 1.10 -0.08 3.64 1.00 2009 2012
## y1[29] -0.73 -0.71 1.15 1.12 -2.66 1.10 1.00 1811 1859
## y1[30] -2.34 -2.33 1.16 1.19 -4.25 -0.48 1.00 1546 1538
## y1[31] 0.97 0.97 1.12 1.09 -0.89 2.78 1.00 2044 1917
## y1[32] 0.08 0.10 1.14 1.13 -1.81 1.89 1.00 2023 1736
## y1[33] 0.94 0.95 1.14 1.14 -1.00 2.79 1.00 1896 1840
## y1[34] -1.38 -1.37 1.15 1.17 -3.24 0.50 1.00 2059 1821
## y1[35] 3.19 3.17 1.17 1.14 1.26 5.16 1.00 1986 1894
## y1[36] 2.94 2.91 1.17 1.13 1.04 4.88 1.00 1918 1927
## y1[37] -1.18 -1.19 1.17 1.18 -3.08 0.78 1.00 2214 1970
## y1[38] 0.46 0.47 1.13 1.12 -1.42 2.32 1.00 1873 2055
## y1[39] 0.48 0.45 1.12 1.09 -1.34 2.37 1.00 1933 1926
## y1[40] 0.50 0.48 1.15 1.09 -1.37 2.43 1.00 1696 2056
## y1[41] 1.58 1.60 1.13 1.17 -0.31 3.34 1.00 1970 1968
## y1[42] 2.78 2.78 1.13 1.14 0.93 4.59 1.00 1886 1776
## y1[43] -1.86 -1.86 1.20 1.19 -3.77 0.06 1.00 1417 1883
## y1[44] 0.36 0.36 1.16 1.15 -1.51 2.33 1.00 1544 1885
## y1[45] -0.92 -0.91 1.14 1.14 -2.73 0.92 1.00 2105 1930
## y1[46] 1.16 1.17 1.13 1.12 -0.70 3.00 1.00 2025 1959
## y1[47] -2.71 -2.74 1.20 1.16 -4.65 -0.66 1.00 1632 2008
## y1[48] 2.91 2.91 1.23 1.22 0.86 4.96 1.00 1766 1793
## y1[49] 3.86 3.87 1.18 1.15 1.86 5.70 1.00 1639 1841
## y1[50] 1.24 1.27 1.18 1.14 -0.69 3.08 1.00 1673 1953
## m1[1] -1.39 -1.40 0.31 0.30 -1.89 -0.88 1.00 1698 1677
## m1[2] -1.53 -1.54 0.31 0.31 -2.04 -1.00 1.00 1724 1664
## m1[3] 0.55 0.55 0.37 0.38 -0.06 1.14 1.00 879 907
## m1[4] 1.11 1.11 0.28 0.28 0.63 1.56 1.00 1659 1369
## m1[5] 1.52 1.52 0.29 0.28 1.05 1.98 1.00 1434 1504
## m1[6] -0.52 -0.52 0.27 0.27 -0.96 -0.06 1.00 1551 1594
## m1[7] 1.01 1.01 0.28 0.28 0.54 1.46 1.00 1651 1295
## m1[8] 1.32 1.32 0.36 0.37 0.73 1.92 1.00 862 891
## m1[9] 3.19 3.20 0.34 0.33 2.63 3.73 1.00 2021 1591
## m1[10] 2.57 2.58 0.31 0.30 2.06 3.07 1.00 1912 1619
## m1[11] 3.47 3.48 0.36 0.34 2.88 4.03 1.00 2048 1671
## m1[12] 3.40 3.40 0.39 0.38 2.75 4.01 1.00 1595 1614
## m1[13] -0.82 -0.83 0.28 0.28 -1.28 -0.36 1.00 1598 1587
## m1[14] 2.22 2.22 0.42 0.42 1.51 2.91 1.00 666 1244
## m1[15] 3.52 3.52 0.42 0.41 2.85 4.22 1.00 1059 1230
## m1[16] -1.51 -1.52 0.31 0.31 -2.02 -0.99 1.00 1721 1664
## m1[17] -0.69 -0.70 0.28 0.28 -1.15 -0.24 1.00 1578 1520
## m1[18] 1.68 1.68 0.37 0.38 1.09 2.28 1.00 867 833
## m1[19] 0.40 0.40 0.37 0.39 -0.22 1.00 1.00 887 1067
## m1[20] 0.94 0.95 0.27 0.27 0.49 1.39 1.00 1423 1425
## m1[21] 3.08 3.07 0.40 0.40 2.44 3.74 1.00 1007 1200
## m1[22] -1.33 -1.32 0.44 0.46 -2.04 -0.62 1.00 1060 1446
## m1[23] 0.19 0.18 0.37 0.36 -0.40 0.77 1.00 540 828
## m1[24] -2.61 -2.61 0.43 0.43 -3.31 -1.89 1.00 697 1088
## m1[25] -1.17 -1.17 0.30 0.29 -1.65 -0.67 1.00 1660 1583
## m1[26] 2.85 2.86 0.46 0.44 2.09 3.59 1.00 736 1332
## m1[27] 0.29 0.28 0.29 0.29 -0.19 0.77 1.00 1577 1465
## m1[28] 1.77 1.77 0.29 0.28 1.29 2.22 1.00 1756 1541
## m1[29] -0.73 -0.73 0.33 0.32 -1.27 -0.20 1.00 1593 1451
## m1[30] -2.36 -2.36 0.42 0.42 -3.05 -1.66 1.00 667 987
## m1[31] 0.99 0.99 0.28 0.28 0.52 1.45 1.00 1650 1312
## m1[32] 0.08 0.08 0.30 0.29 -0.41 0.56 1.00 1568 1420
## m1[33] 0.96 0.96 0.27 0.27 0.51 1.41 1.00 1421 1425
## m1[34] -1.40 -1.40 0.36 0.35 -1.98 -0.81 1.00 1674 1595
## m1[35] 3.22 3.23 0.34 0.33 2.66 3.76 1.00 2022 1617
## m1[36] 2.95 2.96 0.33 0.31 2.42 3.48 1.00 1990 1493
## m1[37] -1.22 -1.22 0.35 0.34 -1.79 -0.66 1.00 1649 1581
## m1[38] 0.47 0.47 0.26 0.26 0.03 0.91 1.00 1438 1542
## m1[39] 0.50 0.49 0.29 0.28 0.02 0.96 1.00 1595 1444
## m1[40] 0.47 0.47 0.26 0.26 0.03 0.91 1.00 1438 1542
## m1[41] 1.59 1.60 0.29 0.29 1.11 2.06 1.00 1438 1564
## m1[42] 2.76 2.75 0.39 0.39 2.14 3.40 1.00 957 1078
## m1[43] -1.89 -1.89 0.40 0.40 -2.54 -1.22 1.00 618 1010
## m1[44] 0.36 0.35 0.37 0.37 -0.24 0.95 1.00 543 836
## m1[45] -0.89 -0.89 0.29 0.28 -1.35 -0.42 1.00 1610 1535
## m1[46] 1.15 1.15 0.28 0.28 0.67 1.60 1.00 1663 1403
## m1[47] -2.74 -2.74 0.44 0.43 -3.45 -2.01 1.00 714 1101
## m1[48] 2.93 2.94 0.46 0.45 2.16 3.67 1.00 745 1370
## m1[49] 3.85 3.86 0.42 0.40 3.15 4.52 1.00 1631 1613
## m1[50] 1.27 1.26 0.36 0.37 0.67 1.87 1.00 861 891
tb=tibble(xa=runif(n,-2,2),xb=runif(n,-2,2),
ca=sample(c('a','b'),n,replace=T),cb=sample(c('a','b'),n,replace=T),
y=rnorm(n,xa+xb-xa*xb+(ca=='b')*2+(cb=='b')*2-(ca=='b')*(cb=='b'),1))
grid.arrange(qplot(data=tb,xa,y,col=ca),
qplot(data=tb,xa,y,col=cb),
qplot(data=tb,xb,y,col=ca),
qplot(data=tb,xb,y,col=cb))
fn=function(f){
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mle=mdl$optimize(data=data)
mle
mcmc=goMCMC(mdl,data)
seeMCMC(mcmc,exc='m',ch=F)
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
grid.arrange(
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
grid.arrange(
qplot(data=tb,1:n,y,col=ca)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
qplot(data=tb,1:n,y,col=cb)+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray'),
ncol=2
)
}
f0=formula(y~xa+xb+ca+cb)
f1=formula(y~xa+xb+ca*cb)
f2=formula(y~xa*xb+ca*cb)
fn(f0)## Initial log joint probability = -847.296
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -48.2121 0.000112663 0.0012299 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -50.93 -50.60 1.81 1.70 -54.47 -48.60 1.00 870 1288
## b[1] -0.02 -0.01 0.41 0.38 -0.70 0.65 1.00 660 812
## b[2] 1.25 1.25 0.20 0.20 0.92 1.57 1.00 2716 1608
## b[3] 1.09 1.10 0.23 0.24 0.72 1.46 1.00 2081 1523
## b[4] 2.04 2.03 0.48 0.46 1.23 2.83 1.00 733 785
## b[5] 1.50 1.48 0.51 0.50 0.68 2.33 1.01 701 856
## s 1.72 1.71 0.19 0.17 1.44 2.04 1.00 1965 1437
## y1[1] -3.26 -3.22 1.83 1.85 -6.20 -0.30 1.00 1881 1834
## y1[2] 0.65 0.66 1.80 1.72 -2.22 3.72 1.00 2032 1895
## y1[3] 1.79 1.81 1.79 1.73 -1.20 4.69 1.00 1963 1641
## y1[4] -0.81 -0.83 1.87 1.80 -3.94 2.26 1.00 1885 1920
## y1[5] 3.20 3.16 1.80 1.75 0.25 6.21 1.00 2000 1969
## y1[6] 4.87 4.86 1.74 1.67 2.03 7.79 1.00 1963 1996
## y1[7] -2.14 -2.07 1.80 1.81 -5.13 0.82 1.00 1784 1932
## y1[8] 0.97 0.97 1.81 1.76 -1.89 4.00 1.00 1931 1892
## y1[9] 3.02 2.95 1.82 1.79 0.06 5.96 1.00 1899 1970
## y1[10] 0.41 0.40 1.84 1.82 -2.67 3.49 1.00 1959 1954
## y1[11] -2.39 -2.43 1.83 1.79 -5.43 0.67 1.00 1957 1953
## y1[12] 3.71 3.77 1.87 1.83 0.72 6.67 1.00 2056 2019
## y1[13] 5.19 5.20 1.87 1.83 2.11 8.31 1.00 1965 1971
## y1[14] 1.13 1.12 1.74 1.70 -1.64 4.07 1.00 1966 1819
## y1[15] -0.59 -0.60 1.84 1.79 -3.60 2.58 1.00 2039 1812
## y1[16] -2.52 -2.53 1.84 1.75 -5.49 0.52 1.00 1761 1766
## y1[17] 4.12 4.13 1.80 1.72 1.18 7.08 1.00 1521 1629
## y1[18] 0.58 0.60 1.79 1.78 -2.39 3.49 1.00 2013 1972
## y1[19] 1.29 1.31 1.80 1.82 -1.74 4.19 1.00 1944 2015
## y1[20] -0.18 -0.18 1.78 1.79 -3.07 2.86 1.00 1971 1951
## y1[21] -3.08 -3.05 1.76 1.70 -5.93 -0.25 1.00 1840 1573
## y1[22] 2.36 2.35 1.91 1.90 -0.80 5.48 1.00 1894 2014
## y1[23] 5.69 5.67 1.88 1.85 2.59 8.73 1.00 1946 1863
## y1[24] 3.12 3.12 1.84 1.81 0.11 6.14 1.00 1704 1963
## y1[25] -0.08 -0.06 1.77 1.69 -3.04 2.81 1.00 1738 1878
## y1[26] -0.45 -0.44 1.79 1.75 -3.35 2.50 1.00 1827 1717
## y1[27] -0.12 -0.08 1.82 1.77 -3.10 2.82 1.00 1507 1890
## y1[28] 2.71 2.72 1.80 1.76 -0.26 5.69 1.00 1764 1664
## y1[29] 1.89 1.88 1.81 1.81 -1.18 4.80 1.00 1778 1709
## y1[30] 5.94 5.97 1.89 1.83 2.79 9.04 1.00 2017 1966
## y1[31] 2.21 2.19 1.89 1.79 -0.85 5.35 1.00 1904 2047
## y1[32] 3.99 4.00 1.87 1.91 0.96 7.04 1.00 2076 1974
## y1[33] -1.19 -1.20 1.76 1.75 -4.05 1.74 1.00 2158 2093
## y1[34] -1.14 -1.14 1.88 1.87 -4.12 1.93 1.00 2050 1913
## y1[35] 2.30 2.31 1.82 1.83 -0.67 5.26 1.00 1953 1879
## y1[36] 5.96 5.95 1.82 1.79 2.96 8.97 1.00 2173 1880
## y1[37] 2.09 2.11 1.77 1.76 -0.89 4.97 1.00 1989 1944
## y1[38] 5.35 5.34 1.77 1.80 2.49 8.26 1.00 1934 1919
## y1[39] 0.95 0.92 1.89 1.76 -2.08 4.11 1.00 1997 2014
## y1[40] -0.48 -0.47 1.76 1.81 -3.36 2.40 1.00 1602 1733
## y1[41] 0.87 0.81 1.80 1.76 -2.14 3.82 1.00 1582 1891
## y1[42] -0.26 -0.27 1.78 1.72 -3.11 2.67 1.00 2019 1822
## y1[43] -2.14 -2.11 1.78 1.73 -5.22 0.76 1.00 1864 1833
## y1[44] 5.06 5.05 1.77 1.75 2.17 7.93 1.00 1920 1991
## y1[45] 0.30 0.31 1.80 1.74 -2.58 3.34 1.00 2116 1918
## y1[46] -3.69 -3.72 1.82 1.76 -6.73 -0.70 1.00 2000 1687
## y1[47] 2.41 2.43 1.79 1.76 -0.45 5.44 1.00 1526 1832
## y1[48] 0.12 0.07 1.81 1.79 -2.83 3.05 1.00 1762 1787
## y1[49] 4.89 4.91 1.78 1.72 1.94 7.69 1.00 1937 1934
## y1[50] 1.49 1.51 1.87 1.74 -1.64 4.67 1.00 1883 1722
## m1[1] -3.21 -3.22 0.51 0.52 -4.04 -2.33 1.00 1301 1239
## m1[2] 0.68 0.68 0.61 0.60 -0.32 1.72 1.00 1366 1273
## m1[3] 1.84 1.85 0.49 0.49 1.07 2.63 1.00 1385 1536
## m1[4] -0.77 -0.79 0.57 0.57 -1.70 0.19 1.00 1405 1431
## m1[5] 3.29 3.29 0.49 0.48 2.49 4.07 1.00 1267 1419
## m1[6] 4.92 4.91 0.47 0.45 4.17 5.68 1.00 1248 1315
## m1[7] -2.11 -2.11 0.51 0.50 -2.96 -1.26 1.00 1219 1311
## m1[8] 0.92 0.92 0.49 0.49 0.13 1.74 1.00 1167 1335
## m1[9] 2.97 2.97 0.60 0.59 1.95 3.96 1.00 956 1215
## m1[10] 0.41 0.41 0.67 0.65 -0.71 1.49 1.00 1740 1474
## m1[11] -2.39 -2.39 0.53 0.52 -3.27 -1.54 1.00 1181 1391
## m1[12] 3.73 3.72 0.55 0.56 2.84 4.61 1.00 1294 1290
## m1[13] 5.23 5.22 0.53 0.52 4.37 6.10 1.00 1433 1325
## m1[14] 1.09 1.08 0.46 0.46 0.37 1.83 1.00 1164 1267
## m1[15] -0.56 -0.57 0.62 0.61 -1.59 0.46 1.00 1565 1520
## m1[16] -2.55 -2.55 0.50 0.49 -3.37 -1.71 1.00 1134 1297
## m1[17] 4.14 4.14 0.56 0.54 3.21 5.06 1.00 1469 1433
## m1[18] 0.61 0.60 0.44 0.42 -0.11 1.34 1.00 1257 1307
## m1[19] 1.29 1.28 0.43 0.43 0.59 2.00 1.00 1188 1245
## m1[20] -0.22 -0.23 0.47 0.45 -0.98 0.55 1.00 1419 1343
## m1[21] -3.03 -3.03 0.52 0.51 -3.84 -2.15 1.00 1310 1305
## m1[22] 2.39 2.39 0.69 0.68 1.22 3.51 1.00 1352 1332
## m1[23] 5.77 5.78 0.66 0.66 4.67 6.85 1.00 1759 1655
## m1[24] 3.12 3.14 0.64 0.63 2.04 4.16 1.00 938 1287
## m1[25] 0.00 0.00 0.50 0.49 -0.83 0.82 1.00 956 1007
## m1[26] -0.45 -0.47 0.47 0.47 -1.22 0.32 1.00 1417 1340
## m1[27] -0.11 -0.11 0.66 0.67 -1.17 0.98 1.00 1273 1537
## m1[28] 2.67 2.66 0.49 0.49 1.87 3.48 1.00 1108 1402
## m1[29] 1.86 1.86 0.48 0.48 1.09 2.65 1.00 1171 1267
## m1[30] 5.89 5.90 0.71 0.73 4.74 7.04 1.00 1374 1525
## m1[31] 2.21 2.21 0.65 0.66 1.16 3.24 1.00 1875 1636
## m1[32] 4.04 4.04 0.51 0.52 3.20 4.87 1.00 1242 1242
## m1[33] -1.18 -1.20 0.55 0.53 -2.09 -0.28 1.00 1784 1410
## m1[34] -1.13 -1.15 0.55 0.53 -2.01 -0.23 1.00 1554 1305
## m1[35] 2.35 2.34 0.58 0.56 1.38 3.30 1.00 1202 1183
## m1[36] 5.97 5.95 0.52 0.50 5.12 6.85 1.00 1524 1378
## m1[37] 2.10 2.10 0.45 0.45 1.39 2.83 1.00 1119 1249
## m1[38] 5.30 5.30 0.49 0.47 4.51 6.11 1.00 1338 1343
## m1[39] 0.91 0.90 0.62 0.61 -0.15 1.90 1.00 1649 1431
## m1[40] -0.44 -0.44 0.43 0.40 -1.14 0.25 1.00 708 864
## m1[41] 0.87 0.88 0.58 0.56 -0.08 1.83 1.00 1101 1327
## m1[42] -0.31 -0.32 0.46 0.46 -1.09 0.46 1.00 897 977
## m1[43] -2.05 -2.06 0.44 0.42 -2.78 -1.33 1.00 914 1116
## m1[44] 5.08 5.06 0.48 0.47 4.30 5.89 1.00 1319 1350
## m1[45] 0.27 0.27 0.44 0.44 -0.46 1.00 1.00 1275 1301
## m1[46] -3.66 -3.67 0.55 0.53 -4.55 -2.74 1.00 1460 1278
## m1[47] 2.38 2.38 0.59 0.57 1.38 3.33 1.00 874 1222
## m1[48] 0.11 0.11 0.65 0.64 -0.95 1.21 1.00 1294 1563
## m1[49] 4.84 4.82 0.49 0.49 4.04 5.68 1.00 1314 1327
## m1[50] 1.51 1.50 0.72 0.71 0.31 2.73 1.00 1447 1380
fn(f1)## Initial log joint probability = -116.493
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 27 -45.1234 0.000147958 0.000662463 1 1 31
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -48.59 -48.23 2.09 1.84 -52.59 -45.95 1.01 695 957
## b[1] -0.51 -0.52 0.46 0.48 -1.24 0.21 1.01 655 871
## b[2] 1.29 1.29 0.19 0.19 0.98 1.59 1.00 2970 1618
## b[3] 1.04 1.04 0.22 0.22 0.66 1.39 1.00 2333 1506
## b[4] 3.02 3.02 0.63 0.63 2.02 4.04 1.01 570 1035
## b[5] 2.58 2.57 0.69 0.71 1.50 3.71 1.01 571 761
## b[6] -2.19 -2.16 0.97 0.98 -3.81 -0.60 1.02 495 590
## s 1.65 1.64 0.19 0.18 1.38 1.97 1.00 1559 1317
## y1[1] -3.68 -3.66 1.71 1.71 -6.54 -0.84 1.00 1791 1770
## y1[2] 1.40 1.34 1.76 1.74 -1.41 4.33 1.00 1647 1932
## y1[3] 2.45 2.45 1.77 1.76 -0.42 5.32 1.00 2145 2055
## y1[4] -0.24 -0.23 1.73 1.68 -3.12 2.63 1.00 1890 1925
## y1[5] 3.91 3.89 1.79 1.78 0.94 6.89 1.00 1811 1973
## y1[6] 4.31 4.34 1.75 1.69 1.43 7.18 1.00 1895 2009
## y1[7] -2.46 -2.45 1.75 1.77 -5.39 0.34 1.00 1896 1746
## y1[8] 1.58 1.57 1.75 1.72 -1.27 4.38 1.00 1853 1919
## y1[9] 2.52 2.57 1.76 1.80 -0.36 5.38 1.00 1805 1909
## y1[10] 0.88 0.83 1.80 1.73 -1.99 3.91 1.00 1875 1916
## y1[11] -2.94 -2.92 1.74 1.75 -5.72 -0.07 1.00 1906 1861
## y1[12] 3.22 3.23 1.76 1.75 0.37 6.18 1.00 2011 2083
## y1[13] 4.67 4.65 1.79 1.78 1.81 7.63 1.00 1988 1945
## y1[14] 1.67 1.66 1.73 1.69 -1.21 4.52 1.00 1937 1895
## y1[15] -0.05 -0.07 1.75 1.71 -2.95 2.79 1.00 2034 2011
## y1[16] -3.09 -3.13 1.75 1.70 -5.87 -0.27 1.00 1724 1811
## y1[17] 4.84 4.81 1.77 1.66 1.98 7.83 1.00 1980 1971
## y1[18] 1.09 1.15 1.73 1.70 -1.80 3.84 1.00 2124 1872
## y1[19] 1.79 1.80 1.67 1.67 -0.96 4.47 1.00 2175 2025
## y1[20] 0.34 0.32 1.69 1.66 -2.41 3.07 1.00 1635 1706
## y1[21] -3.48 -3.45 1.71 1.68 -6.33 -0.71 1.00 1830 1922
## y1[22] 2.66 2.65 1.79 1.74 -0.35 5.68 1.00 1918 1929
## y1[23] 6.31 6.35 1.82 1.74 3.29 9.33 1.00 2054 1969
## y1[24] 2.59 2.63 1.81 1.80 -0.45 5.47 1.00 1743 1995
## y1[25] -0.39 -0.40 1.75 1.75 -3.19 2.53 1.00 2099 1892
## y1[26] 0.04 0.07 1.75 1.64 -2.86 2.90 1.00 1886 1827
## y1[27] -0.67 -0.65 1.81 1.75 -3.64 2.31 1.00 1946 1972
## y1[28] 2.06 2.07 1.73 1.66 -0.76 4.90 1.00 1939 1821
## y1[29] 2.45 2.46 1.75 1.74 -0.43 5.37 1.00 1512 1636
## y1[30] 6.35 6.31 1.79 1.72 3.40 9.21 1.00 1990 1807
## y1[31] 2.92 2.90 1.84 1.75 -0.04 6.05 1.00 1918 2059
## y1[32] 3.53 3.48 1.70 1.71 0.85 6.28 1.00 1697 1892
## y1[33] -0.58 -0.53 1.77 1.76 -3.46 2.36 1.00 2018 2015
## y1[34] -0.68 -0.67 1.76 1.77 -3.65 2.21 1.00 2023 1796
## y1[35] 2.74 2.78 1.74 1.71 -0.16 5.63 1.00 1915 1779
## y1[36] 5.38 5.40 1.72 1.69 2.50 8.24 1.00 1954 1845
## y1[37] 2.69 2.66 1.74 1.74 -0.29 5.46 1.00 1835 1762
## y1[38] 4.66 4.67 1.75 1.76 1.80 7.50 1.00 2082 2007
## y1[39] 1.40 1.40 1.81 1.75 -1.59 4.29 1.00 1985 1894
## y1[40] -0.96 -0.99 1.75 1.79 -3.78 1.78 1.00 1815 1735
## y1[41] 0.59 0.60 1.74 1.69 -2.26 3.50 1.00 1916 1971
## y1[42] -0.75 -0.73 1.70 1.62 -3.48 1.99 1.00 1927 1980
## y1[43] -2.52 -2.50 1.72 1.70 -5.38 0.25 1.00 1884 2008
## y1[44] 4.40 4.41 1.73 1.75 1.53 7.25 1.00 2002 1960
## y1[45] 0.78 0.79 1.70 1.70 -1.91 3.60 1.00 1776 1724
## y1[46] -4.27 -4.28 1.75 1.70 -7.08 -1.39 1.00 1906 1824
## y1[47] 1.84 1.85 1.74 1.63 -1.04 4.70 1.00 1848 1819
## y1[48] -0.53 -0.52 1.77 1.73 -3.46 2.39 1.00 1948 1932
## y1[49] 4.12 4.16 1.74 1.77 1.25 6.90 1.00 2070 1786
## y1[50] 1.85 1.80 1.77 1.70 -0.97 4.78 1.00 1872 1827
## m1[1] -3.69 -3.69 0.52 0.52 -4.52 -2.82 1.01 960 1310
## m1[2] 1.39 1.37 0.69 0.69 0.33 2.57 1.01 1134 1222
## m1[3] 2.43 2.43 0.54 0.53 1.55 3.33 1.00 1431 1533
## m1[4] -0.24 -0.25 0.59 0.59 -1.17 0.73 1.00 1305 1350
## m1[5] 3.91 3.90 0.57 0.57 3.00 4.88 1.00 1135 1235
## m1[6] 4.30 4.29 0.54 0.52 3.42 5.16 1.00 1578 1204
## m1[7] -2.51 -2.52 0.53 0.54 -3.38 -1.65 1.00 996 1622
## m1[8] 1.57 1.55 0.57 0.57 0.68 2.50 1.01 1024 1149
## m1[9] 2.52 2.52 0.63 0.64 1.48 3.52 1.01 1075 1063
## m1[10] 0.84 0.82 0.67 0.65 -0.20 1.93 1.00 1721 1487
## m1[11] -2.96 -2.95 0.54 0.53 -3.85 -2.07 1.01 854 1000
## m1[12] 3.22 3.21 0.60 0.60 2.26 4.20 1.00 1834 1507
## m1[13] 4.69 4.69 0.58 0.57 3.74 5.67 1.00 1820 1468
## m1[14] 1.65 1.64 0.51 0.51 0.83 2.49 1.00 1090 1419
## m1[15] -0.08 -0.09 0.62 0.61 -1.04 0.94 1.00 1483 1443
## m1[16] -3.09 -3.09 0.51 0.50 -3.92 -2.26 1.01 823 969
## m1[17] 4.79 4.79 0.64 0.62 3.77 5.88 1.00 1262 1342
## m1[18] 1.15 1.15 0.48 0.47 0.37 1.95 1.00 1351 1590
## m1[19] 1.83 1.84 0.48 0.46 1.06 2.62 1.00 1322 1433
## m1[20] 0.31 0.32 0.50 0.50 -0.51 1.16 1.00 1461 1520
## m1[21] -3.48 -3.48 0.53 0.54 -4.33 -2.59 1.00 990 1515
## m1[22] 2.74 2.74 0.67 0.64 1.65 3.81 1.00 1667 1665
## m1[23] 6.34 6.33 0.71 0.68 5.16 7.49 1.00 1606 1455
## m1[24] 2.60 2.60 0.67 0.68 1.50 3.67 1.01 1034 1055
## m1[25] -0.39 -0.39 0.52 0.53 -1.25 0.47 1.00 937 1344
## m1[26] 0.06 0.06 0.50 0.49 -0.74 0.89 1.00 1504 1572
## m1[27] -0.70 -0.69 0.69 0.66 -1.88 0.41 1.00 1809 1540
## m1[28] 2.01 2.01 0.55 0.56 1.09 2.90 1.00 1531 1405
## m1[29] 2.51 2.49 0.57 0.55 1.64 3.43 1.01 1024 1110
## m1[30] 6.39 6.40 0.71 0.69 5.20 7.48 1.00 1612 1519
## m1[31] 2.88 2.87 0.70 0.68 1.73 4.03 1.00 1710 1796
## m1[32] 3.51 3.50 0.57 0.57 2.58 4.44 1.00 1747 1638
## m1[33] -0.62 -0.61 0.58 0.58 -1.59 0.36 1.00 1632 1417
## m1[34] -0.68 -0.70 0.54 0.52 -1.54 0.25 1.00 1715 1509
## m1[35] 2.75 2.74 0.57 0.55 1.81 3.68 1.00 1492 1553
## m1[36] 5.37 5.38 0.58 0.55 4.44 6.34 1.00 1735 1369
## m1[37] 2.71 2.71 0.53 0.53 1.88 3.59 1.00 1022 1179
## m1[38] 4.71 4.70 0.55 0.53 3.81 5.63 1.00 1653 1398
## m1[39] 1.36 1.34 0.63 0.61 0.37 2.42 1.00 1597 1360
## m1[40] -0.96 -0.97 0.46 0.47 -1.72 -0.22 1.01 657 975
## m1[41] 0.51 0.52 0.60 0.61 -0.48 1.49 1.00 1190 1713
## m1[42] -0.72 -0.73 0.49 0.51 -1.54 0.08 1.00 853 1280
## m1[43] -2.53 -2.53 0.46 0.47 -3.27 -1.76 1.01 745 1093
## m1[44] 4.44 4.43 0.55 0.53 3.54 5.34 1.00 1583 1288
## m1[45] 0.80 0.81 0.48 0.46 0.02 1.59 1.00 1386 1586
## m1[46] -4.16 -4.16 0.55 0.55 -5.04 -3.26 1.01 1024 1226
## m1[47] 1.84 1.85 0.63 0.65 0.82 2.85 1.01 913 918
## m1[48] -0.54 -0.53 0.68 0.67 -1.68 0.55 1.00 1781 1536
## m1[49] 4.18 4.17 0.56 0.54 3.26 5.09 1.00 1569 1325
## m1[50] 1.85 1.84 0.69 0.64 0.73 2.97 1.00 1777 1671
fn(f2)## Initial log joint probability = -1164.03
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 28 -18.0098 7.85082e-05 0.000763031 1 1 35
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -22.53 -22.19 2.16 2.00 -26.56 -19.71 1.00 745 1205
## b[1] 0.01 0.00 0.28 0.28 -0.44 0.48 1.00 673 985
## b[2] 1.24 1.24 0.11 0.11 1.05 1.42 1.00 2311 1595
## b[3] 1.04 1.04 0.13 0.13 0.84 1.26 1.00 2457 1955
## b[4] 2.37 2.38 0.38 0.38 1.74 2.99 1.01 524 723
## b[5] 1.86 1.86 0.40 0.40 1.19 2.52 1.00 450 757
## b[6] -0.91 -0.91 0.10 0.10 -1.07 -0.74 1.00 2037 1465
## b[7] -1.03 -1.05 0.58 0.59 -1.93 -0.04 1.00 398 620
## s 0.97 0.96 0.11 0.11 0.81 1.17 1.00 1958 1703
## y1[1] -4.81 -4.80 1.01 0.97 -6.53 -3.18 1.00 2124 2002
## y1[2] 2.39 2.42 1.03 1.02 0.65 4.08 1.00 1908 1932
## y1[3] 3.16 3.17 1.05 0.99 1.37 4.88 1.00 1772 1627
## y1[4] -0.74 -0.73 1.04 1.08 -2.43 0.98 1.00 1970 1583
## y1[5] 3.40 3.43 1.04 1.02 1.75 5.05 1.00 1872 1792
## y1[6] 4.25 4.28 1.04 1.02 2.54 5.98 1.00 1979 2035
## y1[7] -2.20 -2.22 1.00 1.02 -3.78 -0.56 1.00 1891 2009
## y1[8] 1.62 1.62 1.04 1.04 -0.07 3.35 1.00 1574 1611
## y1[9] 1.80 1.79 1.05 0.99 0.10 3.57 1.00 1779 1732
## y1[10] 2.85 2.84 1.06 1.01 1.10 4.65 1.00 2110 1930
## y1[11] -2.49 -2.48 1.04 1.07 -4.20 -0.83 1.00 1645 1854
## y1[12] 4.69 4.67 1.05 1.05 3.00 6.43 1.00 2053 1820
## y1[13] 5.34 5.34 1.01 1.01 3.67 6.98 1.00 1879 1884
## y1[14] 1.58 1.56 1.03 1.00 -0.15 3.24 1.00 1781 1702
## y1[15] 0.39 0.41 1.01 0.98 -1.27 2.00 1.00 1871 1860
## y1[16] -3.14 -3.12 1.02 0.99 -4.88 -1.52 1.00 1450 1868
## y1[17] 4.15 4.13 1.04 1.00 2.44 5.92 1.00 2023 1759
## y1[18] 0.89 0.89 1.02 0.97 -0.77 2.60 1.00 1805 1777
## y1[19] 1.84 1.85 1.01 0.99 0.17 3.51 1.00 1883 1873
## y1[20] -0.53 -0.51 1.05 1.03 -2.26 1.17 1.00 1585 1487
## y1[21] -4.34 -4.35 1.02 0.99 -5.97 -2.62 1.00 1857 1876
## y1[22] 4.47 4.46 1.06 1.03 2.79 6.22 1.00 2081 1972
## y1[23] 2.98 3.00 1.14 1.16 1.06 4.83 1.00 1584 1643
## y1[24] 1.40 1.41 1.06 1.05 -0.33 3.12 1.00 2020 2016
## y1[25] 1.01 0.98 1.04 1.03 -0.66 2.79 1.00 1715 1867
## y1[26] -1.02 -1.01 1.03 1.02 -2.69 0.65 1.00 1853 1621
## y1[27] -2.49 -2.46 1.09 1.07 -4.28 -0.72 1.00 2038 2010
## y1[28] 2.42 2.38 1.04 0.99 0.74 4.19 1.00 2000 2016
## y1[29] 2.72 2.71 1.02 1.00 1.04 4.43 1.00 1877 1880
## y1[30] 3.78 3.79 1.09 1.06 1.96 5.52 1.00 1761 1787
## y1[31] 5.63 5.63 1.08 1.07 3.87 7.44 1.00 1839 1854
## y1[32] 4.57 4.59 1.02 0.98 2.89 6.21 1.00 1938 1798
## y1[33] -2.24 -2.26 1.07 1.05 -4.04 -0.50 1.00 1748 1947
## y1[34] -2.00 -1.99 1.01 1.02 -3.63 -0.34 1.00 1719 1729
## y1[35] 3.42 3.43 1.03 0.98 1.75 5.10 1.00 1901 1965
## y1[36] 4.78 4.81 1.03 1.05 3.05 6.47 1.00 1904 1890
## y1[37] 2.52 2.51 1.02 1.01 0.86 4.23 1.00 1639 1916
## y1[38] 4.61 4.61 1.03 0.99 2.97 6.26 1.00 1897 1972
## y1[39] 2.91 2.92 1.03 1.02 1.22 4.63 1.00 2131 1907
## y1[40] -0.36 -0.35 0.97 0.96 -1.93 1.24 1.00 1533 1603
## y1[41] 2.60 2.61 1.05 1.07 0.89 4.34 1.00 1719 1936
## y1[42] 0.33 0.33 1.02 1.01 -1.33 1.96 1.00 1833 1853
## y1[43] -2.64 -2.63 0.99 0.96 -4.31 -1.07 1.00 2122 1650
## y1[44] 4.31 4.30 1.04 1.01 2.64 6.02 1.00 1864 1996
## y1[45] 0.28 0.28 1.02 1.00 -1.35 1.98 1.00 1836 1839
## y1[46] -5.64 -5.66 1.06 1.04 -7.37 -3.91 1.00 1976 2054
## y1[47] 1.42 1.40 1.03 1.04 -0.20 3.15 1.00 1385 1805
## y1[48] -1.76 -1.74 1.06 1.01 -3.46 -0.04 1.00 1781 1834
## y1[49] 4.29 4.29 1.02 1.01 2.68 5.99 1.00 1884 1773
## y1[50] 4.29 4.32 1.11 1.12 2.50 6.07 1.00 1951 1896
## m1[1] -4.79 -4.80 0.33 0.32 -5.33 -4.24 1.00 1428 1627
## m1[2] 2.37 2.37 0.39 0.40 1.73 3.03 1.00 966 1626
## m1[3] 3.15 3.15 0.32 0.32 2.64 3.67 1.00 1091 1566
## m1[4] -0.73 -0.72 0.34 0.33 -1.29 -0.17 1.00 991 1296
## m1[5] 3.40 3.41 0.32 0.31 2.83 3.92 1.00 721 1073
## m1[6] 4.28 4.29 0.31 0.30 3.78 4.79 1.00 1428 1579
## m1[7] -2.19 -2.20 0.31 0.30 -2.69 -1.67 1.00 1161 1240
## m1[8] 1.61 1.61 0.31 0.30 1.10 2.13 1.00 707 1300
## m1[9] 1.82 1.82 0.39 0.37 1.18 2.43 1.00 1080 1164
## m1[10] 2.85 2.85 0.45 0.45 2.11 3.59 1.00 2124 1714
## m1[11] -2.51 -2.52 0.33 0.33 -3.04 -1.97 1.00 810 1101
## m1[12] 4.74 4.75 0.36 0.36 4.14 5.34 1.00 1624 1569
## m1[13] 5.34 5.35 0.32 0.32 4.80 5.88 1.00 1632 1584
## m1[14] 1.54 1.54 0.28 0.28 1.06 1.99 1.00 786 1324
## m1[15] 0.40 0.40 0.36 0.36 -0.19 1.00 1.00 1417 1427
## m1[16] -3.13 -3.14 0.31 0.30 -3.62 -2.61 1.00 856 1284
## m1[17] 4.11 4.12 0.37 0.34 3.47 4.70 1.00 800 1205
## m1[18] 0.84 0.84 0.29 0.28 0.35 1.31 1.01 887 1386
## m1[19] 1.85 1.85 0.28 0.27 1.38 2.32 1.01 884 1343
## m1[20] -0.55 -0.54 0.32 0.31 -1.08 -0.04 1.01 949 1418
## m1[21] -4.36 -4.37 0.32 0.32 -4.88 -3.82 1.00 1440 1543
## m1[22] 4.48 4.48 0.45 0.45 3.76 5.22 1.00 1937 1585
## m1[23] 2.99 3.01 0.56 0.55 2.05 3.88 1.00 1145 1304
## m1[24] 1.42 1.43 0.42 0.41 0.73 2.10 1.00 1103 1272
## m1[25] 1.03 1.03 0.34 0.34 0.48 1.59 1.00 1008 1297
## m1[26] -1.03 -1.03 0.33 0.32 -1.57 -0.50 1.01 991 1507
## m1[27] -2.50 -2.50 0.46 0.45 -3.27 -1.75 1.00 1903 1451
## m1[28] 2.41 2.42 0.33 0.31 1.85 2.93 1.00 1269 1400
## m1[29] 2.70 2.70 0.31 0.30 2.18 3.22 1.00 726 1285
## m1[30] 3.78 3.79 0.52 0.50 2.90 4.63 1.01 1186 1448
## m1[31] 5.60 5.60 0.49 0.50 4.82 6.40 1.00 1970 1856
## m1[32] 4.58 4.59 0.33 0.32 4.04 5.11 1.00 1541 1491
## m1[33] -2.20 -2.21 0.39 0.38 -2.84 -1.55 1.00 1096 1495
## m1[34] -1.99 -1.99 0.36 0.37 -2.57 -1.42 1.01 1221 1400
## m1[35] 3.43 3.42 0.35 0.36 2.86 4.01 1.01 1476 1302
## m1[36] 4.79 4.79 0.34 0.34 4.25 5.36 1.00 1669 1558
## m1[37] 2.53 2.54 0.29 0.28 2.03 3.00 1.00 678 1053
## m1[38] 4.61 4.60 0.31 0.31 4.10 5.12 1.00 1527 1506
## m1[39] 2.92 2.92 0.40 0.40 2.27 3.57 1.00 1942 1658
## m1[40] -0.35 -0.36 0.29 0.28 -0.82 0.12 1.00 663 941
## m1[41] 2.61 2.61 0.42 0.42 1.94 3.29 1.00 1202 1428
## m1[42] 0.38 0.37 0.31 0.32 -0.12 0.89 1.00 922 1192
## m1[43] -2.65 -2.66 0.28 0.28 -3.09 -2.19 1.00 846 1327
## m1[44] 4.32 4.33 0.32 0.31 3.81 4.85 1.00 1443 1563
## m1[45] 0.25 0.25 0.29 0.29 -0.25 0.72 1.01 901 1420
## m1[46] -5.64 -5.65 0.36 0.36 -6.23 -5.03 1.00 1631 1545
## m1[47] 1.45 1.45 0.37 0.36 0.85 2.06 1.00 920 1206
## m1[48] -1.76 -1.76 0.43 0.42 -2.48 -1.04 1.00 1770 1404
## m1[49] 4.26 4.26 0.33 0.32 3.73 4.79 1.00 1401 1480
## m1[50] 4.31 4.30 0.51 0.52 3.48 5.13 1.00 2083 1730
generalized linear regression
log normal regression
objective variable [0,Infinity)
# of samples is N,
log mi=sum(bj*xji),j=[0,K],i=[1,N]
log yi~N(mi,s) ex5-2.stan
log normal regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> s;
}
model{
vector[N] m=X*b;
y~lognormal(m,s);
}
generated quantities{
vector[N] y1;
vector[N] m1=X*b;
for(i in 1:N){
y1[i]=lognormal_rng(m1[i],s);
}
}n=20
tb=tibble(x1=runif(n,0,9),x2=runif(n,0,9),
y=rnorm(n,log(x1+x2),1))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex5-2.stan')
mle=mdl$optimize(data=data) # it sometimes occur error and stop process## Initial log joint probability = -1103.3
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 32 -11.1287 4.18458e-05 0.0013495 0.9083 0.9083 46
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.13
## b[1] 0.81
## b[2] 0.08
## b[3] -0.07
## s 0.42
## y1[1] 1.88
## y1[2] 6.51
## y1[3] 5.03
## y1[4] 3.10
## y1[5] 3.86
## y1[6] 2.22
## y1[7] 3.00
## y1[8] 4.03
## y1[9] 2.58
## y1[10] 2.96
## y1[11] 6.59
## y1[12] 1.35
## y1[13] 2.04
## y1[14] 2.06
## y1[15] 2.73
## y1[16] 3.81
## y1[17] 3.13
## y1[18] 1.81
## y1[19] 0.80
## y1[20] 6.89
## m1[1] 0.87
## m1[2] 1.35
## m1[3] 1.01
## m1[4] 0.72
## m1[5] 0.90
## m1[6] 0.72
## m1[7] 1.00
## m1[8] 1.15
## m1[9] 0.88
## m1[10] 0.78
## m1[11] 1.31
## m1[12] 0.38
## m1[13] 1.28
## m1[14] 0.58
## m1[15] 0.59
## m1[16] 1.11
## m1[17] 1.19
## m1[18] 0.41
## m1[19] 0.21
## m1[20] 1.33mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='m',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -14.30 -13.86 1.74 1.39 -17.72 -12.40 1.01 386 567
## b[1] 0.81 0.81 0.34 0.34 0.23 1.38 1.02 283 360
## b[2] 0.08 0.08 0.04 0.04 0.01 0.15 1.01 466 516
## b[3] -0.07 -0.07 0.04 0.04 -0.14 0.00 1.01 647 614
## s 0.50 0.49 0.10 0.09 0.37 0.71 1.01 486 648
## y1[1] 2.73 2.38 1.68 1.21 0.94 5.69 1.00 2048 1881
## y1[2] 4.54 3.86 2.82 1.98 1.59 9.75 1.00 1774 1890
## y1[3] 3.12 2.74 1.79 1.38 1.12 6.33 1.00 1756 1618
## y1[4] 2.38 2.06 1.61 1.03 0.83 4.98 1.00 2039 1864
## y1[5] 2.83 2.41 1.80 1.20 0.99 5.99 1.00 1883 1935
## y1[6] 2.38 2.10 1.41 1.06 0.82 4.94 1.00 1484 1868
## y1[7] 3.10 2.71 1.79 1.33 1.15 6.17 1.00 1860 1414
## y1[8] 3.60 3.11 2.19 1.56 1.28 7.42 1.00 1984 1666
## y1[9] 2.87 2.45 1.87 1.28 0.95 6.15 1.01 773 844
## y1[10] 2.56 2.15 1.66 1.07 0.89 5.50 1.00 1818 1686
## y1[11] 4.25 3.69 2.50 1.84 1.54 8.59 1.00 1788 1845
## y1[12] 1.71 1.45 1.09 0.77 0.56 3.62 1.00 1409 1283
## y1[13] 4.12 3.60 2.34 1.79 1.46 8.59 1.00 1750 1605
## y1[14] 2.04 1.79 1.14 0.93 0.70 4.16 1.00 1909 1891
## y1[15] 2.05 1.77 1.28 0.90 0.72 4.26 1.00 1785 1654
## y1[16] 3.51 3.05 2.05 1.51 1.24 7.14 1.00 1920 1739
## y1[17] 3.72 3.25 2.16 1.62 1.31 7.45 1.00 1927 1930
## y1[18] 1.76 1.51 1.08 0.76 0.64 3.77 1.00 1892 1880
## y1[19] 1.46 1.25 0.91 0.65 0.48 3.10 1.00 1573 1992
## y1[20] 4.38 3.80 2.91 1.92 1.52 8.94 1.00 1999 2001
## m1[1] 0.86 0.86 0.18 0.17 0.58 1.13 1.00 1075 1204
## m1[2] 1.34 1.35 0.23 0.21 0.96 1.71 1.01 805 971
## m1[3] 1.00 1.00 0.19 0.19 0.69 1.30 1.01 859 1132
## m1[4] 0.71 0.71 0.20 0.18 0.39 1.02 1.00 1100 1234
## m1[5] 0.89 0.89 0.17 0.16 0.62 1.15 1.00 1189 1292
## m1[6] 0.72 0.72 0.18 0.18 0.42 1.01 1.02 336 641
## m1[7] 0.99 1.00 0.13 0.13 0.77 1.20 1.00 2079 1304
## m1[8] 1.14 1.15 0.16 0.15 0.87 1.40 1.00 1679 1030
## m1[9] 0.88 0.87 0.24 0.24 0.48 1.27 1.02 288 298
## m1[10] 0.77 0.76 0.21 0.21 0.42 1.10 1.01 809 1127
## m1[11] 1.31 1.31 0.21 0.20 0.95 1.65 1.01 848 888
## m1[12] 0.37 0.38 0.25 0.24 -0.05 0.76 1.01 516 676
## m1[13] 1.28 1.28 0.18 0.17 0.99 1.56 1.00 1694 1417
## m1[14] 0.57 0.58 0.17 0.17 0.29 0.85 1.01 580 987
## m1[15] 0.59 0.59 0.17 0.16 0.31 0.86 1.01 584 935
## m1[16] 1.10 1.10 0.18 0.17 0.82 1.39 1.00 1176 1155
## m1[17] 1.18 1.18 0.17 0.16 0.90 1.45 1.00 929 1034
## m1[18] 0.41 0.41 0.21 0.19 0.07 0.73 1.00 1844 1032
## m1[19] 0.21 0.22 0.27 0.24 -0.24 0.62 1.00 1187 966
## m1[20] 1.32 1.33 0.21 0.19 0.98 1.66 1.00 1036 941
y1=mcmc$draws('y1')
smy=summary(y1)
tb$y1=smy$median
tb$l5=smy$q5
tb$u5=smy$q95
par(mfrow=c(1,1))
plot(tb$y,tb$y1,xlab='obs.',ylab='prd.')
abline(0,1)
qplot(data=tb,y,y1,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1)
par(mfrow=c(1,2))
hist(tb$y-tb$y1,xlab='obs.-prd.',main='residual')
density(tb$y-tb$y1) |> plot(xlab='obs.-prd.',main='residual')
tb=arrange(tb,y1)
par(mfrow=c(1,1))
ylim=c(min(tb$l5),max(tb$u5))
plot(tb$y,ylim=ylim,ylab='y',col='red')
par(new=T)
plot(tb$y1,ylim=ylim,ylab='',col='black')
par(new=T)
plot(tb$l5,ylim=ylim,ylab='',type='l',col='gray')
par(new=T)
plot(tb$u5,ylim=ylim,ylab='',type='l',col='gray')
qplot(1:n,y,data=tb,col=I('red'))+
geom_point(aes(x=1:n,y=y1),tb,col='black')+
geom_line(aes(x=1:n,y=l5),tb,col='gray')+
geom_line(aes(x=1:n,y=u5),tb,col='gray')
poisson regression
objective variable [0,Infinity), integer. variance of error is near to mean
(normal linear regression, correlation is near to 1,-1, variance of error is constant)
# of samples is N,
log li=sum(bj*xji),j=[0,K],i=[1,N]
yi~Po(li)
or
li=sum(bj*xji),j=[0,k]
yi~Po(exp li)
when xj -> xj +1, y -> y* exp bj ex6-1.stan
poisson regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] l=X*b;
y~poisson_log(l);
}
generated quantities{
array[N] int y1;
vector[N] l1=X*b;
for(i in 1:N){
y1[i]=poisson_log_rng(l1[i]);
}
}n=30
tb=tibble(x=runif(n,-1,1),c=sample(c('a','b'),n,replace=T),
y=rpois(n,exp(x+(c=='b')*0.5)))
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='poisson')##
## Call: glm(formula = f, family = "poisson", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.314 1.136 0.122
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 57.1
## Residual Deviance: 37.2 AIC: 104X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-1.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -23.2474
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 12 -9.7876 9.45541e-05 0.000384669 1 1 14
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -9.79
## b[1] 0.31
## b[2] 1.14
## b[3] 0.12
## y1[1] 1.00
## y1[2] 3.00
## y1[3] 1.00
## y1[4] 2.00
## y1[5] 12.00
## y1[6] 0.00
## y1[7] 3.00
## y1[8] 1.00
## y1[9] 3.00
## y1[10] 2.00
## y1[11] 1.00
## y1[12] 3.00
## y1[13] 0.00
## y1[14] 2.00
## y1[15] 3.00
## y1[16] 0.00
## y1[17] 0.00
## y1[18] 6.00
## y1[19] 0.00
## y1[20] 2.00
## y1[21] 1.00
## y1[22] 1.00
## y1[23] 3.00
## y1[24] 2.00
## y1[25] 1.00
## y1[26] 3.00
## y1[27] 4.00
## y1[28] 1.00
## y1[29] 2.00
## y1[30] 1.00
## l1[1] 0.81
## l1[2] 0.87
## l1[3] 0.75
## l1[4] 0.60
## l1[5] 1.41
## l1[6] 1.04
## l1[7] 0.75
## l1[8] 0.30
## l1[9] 0.38
## l1[10] -0.54
## l1[11] 1.20
## l1[12] 1.45
## l1[13] 0.38
## l1[14] -0.35
## l1[15] -0.05
## l1[16] -0.65
## l1[17] 0.55
## l1[18] 1.42
## l1[19] -0.17
## l1[20] 0.58
## l1[21] -0.13
## l1[22] -0.36
## l1[23] 0.81
## l1[24] 1.00
## l1[25] 0.00
## l1[26] 0.93
## l1[27] 1.44
## l1[28] 0.48
## l1[29] -0.13
## l1[30] -0.69mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='l1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -11.28 -10.98 1.19 0.99 -13.62 -9.97 1.00 855 1191
## b[1] 0.28 0.28 0.23 0.23 -0.12 0.65 1.00 1095 1155
## b[2] 1.14 1.13 0.27 0.26 0.71 1.58 1.00 1183 1108
## b[3] 0.14 0.13 0.27 0.27 -0.30 0.58 1.00 988 1021
## y1[1] 2.22 2.00 1.51 1.48 0.00 5.00 1.00 1810 1893
## y1[2] 2.41 2.00 1.64 1.48 0.00 5.00 1.00 1940 1932
## y1[3] 2.07 2.00 1.47 1.48 0.00 5.00 1.00 2094 1800
## y1[4] 1.79 2.00 1.36 1.48 0.00 4.00 1.00 1902 1875
## y1[5] 4.15 4.00 2.19 2.97 1.00 8.00 1.00 2026 1944
## y1[6] 2.81 3.00 1.76 1.48 0.00 6.00 1.00 2125 2023
## y1[7] 2.12 2.00 1.55 1.48 0.00 5.00 1.00 1724 1854
## y1[8] 1.31 1.00 1.21 1.48 0.00 4.00 1.00 1872 1959
## y1[9] 1.41 1.00 1.21 1.48 0.00 4.00 1.00 2076 2045
## y1[10] 0.57 0.00 0.79 0.00 0.00 2.00 1.00 2004 1884
## y1[11] 3.29 3.00 2.03 1.48 0.00 7.00 1.00 1853 1677
## y1[12] 4.26 4.00 2.26 1.48 1.00 8.00 1.00 2048 2016
## y1[13] 1.50 1.00 1.25 1.48 0.00 4.00 1.00 2018 2097
## y1[14] 0.74 1.00 0.88 1.48 0.00 2.00 1.00 1950 1999
## y1[15] 0.96 1.00 1.03 1.48 0.00 3.00 1.00 1357 1745
## y1[16] 0.53 0.00 0.77 0.00 0.00 2.00 1.00 1982 1767
## y1[17] 1.70 2.00 1.35 1.48 0.00 4.00 1.00 1951 1777
## y1[18] 4.22 4.00 2.35 2.97 1.00 8.00 1.00 1889 1894
## y1[19] 0.89 1.00 0.96 1.48 0.00 3.00 1.00 1881 1913
## y1[20] 1.76 2.00 1.34 1.48 0.00 4.00 1.00 2044 2029
## y1[21] 0.87 1.00 0.95 1.48 0.00 3.00 1.00 1791 1921
## y1[22] 0.70 0.00 0.86 0.00 0.00 2.00 1.00 2008 1920
## y1[23] 2.17 2.00 1.51 1.48 0.00 5.00 1.00 2120 1946
## y1[24] 2.74 3.00 1.66 1.48 0.00 6.00 1.00 1951 1959
## y1[25] 0.99 1.00 1.05 1.48 0.00 3.00 1.00 1891 2086
## y1[26] 2.52 2.00 1.67 1.48 0.00 6.00 1.00 1988 1833
## y1[27] 4.22 4.00 2.29 1.48 1.00 8.00 1.00 1848 1702
## y1[28] 1.65 1.00 1.34 1.48 0.00 4.00 1.00 1837 1976
## y1[29] 0.90 1.00 0.97 1.48 0.00 3.00 1.00 2218 1776
## y1[30] 0.49 0.00 0.76 0.00 0.00 2.00 1.00 1892 1869
## l1[1] 0.79 0.80 0.17 0.17 0.50 1.07 1.00 2073 1541
## l1[2] 0.83 0.84 0.22 0.22 0.46 1.17 1.00 1280 1378
## l1[3] 0.73 0.74 0.17 0.17 0.44 1.01 1.00 2055 1486
## l1[4] 0.58 0.58 0.18 0.18 0.27 0.87 1.00 1941 1584
## l1[5] 1.39 1.39 0.21 0.21 1.03 1.73 1.00 1620 1442
## l1[6] 1.02 1.02 0.18 0.17 0.72 1.30 1.00 1999 1432
## l1[7] 0.71 0.72 0.21 0.22 0.35 1.04 1.00 1236 1359
## l1[8] 0.26 0.27 0.24 0.23 -0.15 0.63 1.00 1091 1128
## l1[9] 0.35 0.36 0.21 0.21 0.00 0.69 1.00 1729 1658
## l1[10] -0.57 -0.55 0.38 0.37 -1.20 0.03 1.00 1300 1177
## l1[11] 1.17 1.18 0.24 0.24 0.76 1.55 1.00 1382 1287
## l1[12] 1.43 1.43 0.22 0.22 1.07 1.78 1.00 1587 1380
## l1[13] 0.36 0.36 0.21 0.21 0.01 0.69 1.00 1733 1658
## l1[14] -0.39 -0.38 0.33 0.33 -0.95 0.15 1.00 1053 1035
## l1[15] -0.09 -0.08 0.28 0.27 -0.57 0.36 1.00 1056 1030
## l1[16] -0.69 -0.67 0.39 0.39 -1.33 -0.06 1.00 1057 1038
## l1[17] 0.52 0.53 0.19 0.19 0.21 0.83 1.00 1895 1659
## l1[18] 1.38 1.40 0.27 0.27 0.94 1.82 1.00 1412 1259
## l1[19] -0.19 -0.18 0.30 0.30 -0.69 0.29 1.00 1387 1311
## l1[20] 0.55 0.56 0.19 0.19 0.24 0.85 1.00 1921 1611
## l1[21] -0.18 -0.16 0.29 0.28 -0.68 0.30 1.00 1052 1042
## l1[22] -0.39 -0.37 0.34 0.33 -0.96 0.15 1.00 1330 1258
## l1[23] 0.78 0.79 0.17 0.17 0.50 1.06 1.00 2072 1542
## l1[24] 0.98 0.98 0.17 0.17 0.68 1.26 1.00 2031 1433
## l1[25] -0.04 -0.03 0.27 0.26 -0.51 0.40 1.00 1058 1041
## l1[26] 0.90 0.91 0.22 0.22 0.53 1.24 1.00 1308 1377
## l1[27] 1.42 1.43 0.22 0.22 1.06 1.77 1.00 1592 1380
## l1[28] 0.45 0.45 0.22 0.22 0.07 0.78 1.00 1139 1248
## l1[29] -0.15 -0.14 0.29 0.29 -0.64 0.32 1.00 1402 1409
## l1[30] -0.73 -0.71 0.39 0.40 -1.38 -0.09 1.00 1057 1051
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4
## 0 3 2 2 0 0
## 1 0 6 2 0 1
## 2 1 1 3 2 0
## 4 0 0 2 1 0
## 5 0 1 0 0 1
## 6 0 0 0 0 1
## 7 0 0 0 0 1cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4
## 0 2 0 1 0 0
## 1 0 5 0 0 0
## 2 0 1 1 0 0
## 4 0 0 1 1 0
## 5 0 0 0 0 1
## 6 0 0 0 0 0
## 7 0 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4
## 0 1 2 1 0 0
## 1 0 1 2 0 1
## 2 1 0 2 2 0
## 4 0 0 1 0 0
## 5 0 1 0 0 0
## 6 0 0 0 0 1
## 7 0 0 0 0 1par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 4
## 0 4 2 1 0 0
## 1 4 3 1 0 1
## 2 1 3 2 1 0
## 4 0 0 2 1 0
## 5 1 0 0 0 1
## 6 0 0 0 0 1
## 7 0 0 0 1 0cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 4
## 0 2 0 1 0 0
## 1 4 1 0 0 0
## 2 0 1 1 0 0
## 4 0 0 1 1 0
## 5 0 0 0 0 1
## 6 0 0 0 0 0
## 7 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 4
## 0 2 2 0 0 0
## 1 0 2 1 0 1
## 2 1 2 1 1 0
## 4 0 0 1 0 0
## 5 1 0 0 0 0
## 6 0 0 0 0 1
## 7 0 0 0 1 0par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
logistic regression
# of samples is N,
objective variable 0/1 binary
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~Ber(pi), 0/1 binary
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-2.stan
logistic regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~bernoulli_logit(z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=bernoulli_rng(inv_logit(z1[i]));
}
}n=30
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,1,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
glm(f,tb,family='binomial') # it can caluculte when all trials are once##
## Call: glm(formula = f, family = "binomial", data = tb)
##
## Coefficients:
## (Intercept) x cb
## 0.33 1.03 1.46
##
## Degrees of Freedom: 29 Total (i.e. Null); 27 Residual
## Null Deviance: 32.6
## Residual Deviance: 30.1 AIC: 36.1X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-2.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -33.4797
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 10 -15.0554 0.000209731 0.000287005 0.7124 0.7124 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -15.06
## b[1] 0.33
## b[2] 1.03
## b[3] 1.46
## y1[1] 0.00
## y1[2] 1.00
## y1[3] 0.00
## y1[4] 1.00
## y1[5] 1.00
## y1[6] 0.00
## y1[7] 1.00
## y1[8] 0.00
## y1[9] 1.00
## y1[10] 0.00
## y1[11] 1.00
## y1[12] 0.00
## y1[13] 1.00
## y1[14] 0.00
## y1[15] 1.00
## y1[16] 1.00
## y1[17] 1.00
## y1[18] 0.00
## y1[19] 1.00
## y1[20] 0.00
## y1[21] 1.00
## y1[22] 0.00
## y1[23] 1.00
## y1[24] 1.00
## y1[25] 1.00
## y1[26] 1.00
## y1[27] 1.00
## y1[28] 0.00
## y1[29] 1.00
## y1[30] 0.00
## z1[1] 1.28
## z1[2] 0.57
## z1[3] 1.16
## z1[4] 2.55
## z1[5] 2.33
## z1[6] 1.72
## z1[7] 0.91
## z1[8] 0.51
## z1[9] 0.68
## z1[10] 1.29
## z1[11] 1.39
## z1[12] 0.63
## z1[13] 2.36
## z1[14] 1.08
## z1[15] 2.07
## z1[16] 0.80
## z1[17] 1.17
## z1[18] 1.18
## z1[19] 2.67
## z1[20] 2.76
## z1[21] 1.24
## z1[22] 1.27
## z1[23] 2.20
## z1[24] 1.08
## z1[25] 1.01
## z1[26] 1.22
## z1[27] 0.52
## z1[28] 0.68
## z1[29] 1.97
## z1[30] -0.54mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -16.74 -16.36 1.39 1.11 -19.50 -15.24 1.00 678 776
## b[1] 0.28 0.26 0.77 0.73 -0.96 1.55 1.00 873 916
## b[2] 1.29 1.24 1.11 1.07 -0.41 3.23 1.01 786 551
## b[3] 1.86 1.78 1.27 1.17 -0.13 4.04 1.00 697 628
## y1[1] 0.79 1.00 0.41 0.00 0.00 1.00 1.00 1653 NA
## y1[2] 0.63 1.00 0.48 0.00 0.00 1.00 1.00 1912 NA
## y1[3] 0.77 1.00 0.42 0.00 0.00 1.00 1.00 2112 NA
## y1[4] 0.92 1.00 0.27 0.00 0.00 1.00 1.00 1780 NA
## y1[5] 0.91 1.00 0.28 0.00 0.00 1.00 1.00 1790 NA
## y1[6] 0.84 1.00 0.36 0.00 0.00 1.00 1.00 1941 NA
## y1[7] 0.74 1.00 0.44 0.00 0.00 1.00 1.00 2054 NA
## y1[8] 0.60 1.00 0.49 0.00 0.00 1.00 1.00 1793 NA
## y1[9] 0.65 1.00 0.48 0.00 0.00 1.00 1.00 1996 NA
## y1[10] 0.79 1.00 0.41 0.00 0.00 1.00 1.00 2002 NA
## y1[11] 0.81 1.00 0.40 0.00 0.00 1.00 1.00 1942 NA
## y1[12] 0.66 1.00 0.47 0.00 0.00 1.00 1.00 1799 NA
## y1[13] 0.92 1.00 0.28 0.00 0.00 1.00 1.00 1630 NA
## y1[14] 0.75 1.00 0.43 0.00 0.00 1.00 1.00 1824 NA
## y1[15] 0.90 1.00 0.30 0.00 0.00 1.00 1.00 1837 NA
## y1[16] 0.68 1.00 0.47 0.00 0.00 1.00 1.00 1926 NA
## y1[17] 0.76 1.00 0.42 0.00 0.00 1.00 1.00 2035 NA
## y1[18] 0.78 1.00 0.42 0.00 0.00 1.00 1.00 1997 NA
## y1[19] 0.92 1.00 0.27 0.00 0.00 1.00 1.00 1999 NA
## y1[20] 0.93 1.00 0.26 0.00 0.00 1.00 1.00 1484 NA
## y1[21] 0.79 1.00 0.41 0.00 0.00 1.00 1.00 1874 NA
## y1[22] 0.80 1.00 0.40 0.00 0.00 1.00 1.00 1998 NA
## y1[23] 0.91 1.00 0.28 0.00 0.00 1.00 1.00 1790 NA
## y1[24] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 1738 NA
## y1[25] 0.74 1.00 0.44 0.00 0.00 1.00 1.00 1789 NA
## y1[26] 0.76 1.00 0.43 0.00 0.00 1.00 1.00 2004 NA
## y1[27] 0.60 1.00 0.49 0.00 0.00 1.00 1.00 1878 NA
## y1[28] 0.67 1.00 0.47 0.00 0.00 1.00 1.00 1897 NA
## y1[29] 0.89 1.00 0.32 0.00 0.00 1.00 1.00 1697 NA
## y1[30] 0.36 0.00 0.48 0.00 0.00 1.00 1.00 1518 NA
## z1[1] 1.47 1.42 0.80 0.79 0.24 2.86 1.00 1555 1475
## z1[2] 0.58 0.57 0.64 0.63 -0.43 1.65 1.00 1171 1161
## z1[3] 1.35 1.28 1.06 0.96 -0.32 3.17 1.00 1671 1286
## z1[4] 3.08 2.97 1.33 1.29 1.17 5.49 1.00 696 728
## z1[5] 2.81 2.70 1.17 1.12 1.11 4.90 1.00 742 746
## z1[6] 2.05 1.98 0.92 0.83 0.69 3.66 1.00 1216 1259
## z1[7] 1.01 0.97 0.61 0.60 0.06 2.08 1.00 1830 1483
## z1[8] 0.50 0.49 0.67 0.65 -0.55 1.63 1.00 1070 1044
## z1[9] 0.71 0.70 0.61 0.60 -0.23 1.73 1.00 1358 1296
## z1[10] 1.48 1.43 0.80 0.80 0.25 2.87 1.00 1544 1474
## z1[11] 1.64 1.56 0.96 0.85 0.20 3.34 1.00 1692 1364
## z1[12] 0.66 0.65 0.62 0.60 -0.31 1.70 1.00 1274 1284
## z1[13] 2.85 2.74 1.19 1.15 1.13 4.97 1.00 733 761
## z1[14] 1.22 1.18 0.68 0.66 0.17 2.42 1.00 1796 1373
## z1[15] 2.49 2.39 1.02 0.95 1.03 4.34 1.00 842 806
## z1[16] 0.86 0.84 0.60 0.59 -0.07 1.87 1.00 1639 1345
## z1[17] 1.33 1.30 0.73 0.72 0.22 2.60 1.00 1688 1483
## z1[18] 1.34 1.31 0.73 0.73 0.22 2.62 1.00 1674 1506
## z1[19] 3.24 3.12 1.43 1.37 1.15 5.80 1.00 683 722
## z1[20] 3.35 3.23 1.51 1.44 1.14 6.05 1.00 674 705
## z1[21] 1.45 1.37 1.01 0.93 -0.13 3.23 1.00 1695 1341
## z1[22] 1.46 1.42 0.79 0.79 0.24 2.84 1.00 1562 1475
## z1[23] 2.65 2.55 1.09 1.03 1.11 4.62 1.00 785 802
## z1[24] 1.24 1.17 1.10 1.01 -0.50 3.14 1.00 1639 1220
## z1[25] 1.13 1.10 0.64 0.62 0.13 2.26 1.00 1847 1485
## z1[26] 1.42 1.34 1.03 0.93 -0.19 3.20 1.00 1689 1304
## z1[27] 0.52 0.51 0.66 0.65 -0.52 1.63 1.00 1098 1090
## z1[28] 0.72 0.71 0.60 0.59 -0.22 1.74 1.00 1374 1314
## z1[29] 2.37 2.26 0.98 0.90 0.96 4.08 1.00 912 941
## z1[30] -0.82 -0.75 1.55 1.48 -3.42 1.64 1.01 718 639
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1
## 0 1 6
## 1 0 23cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1
## 0 1 4
## 1 0 11
##
## , , = b
##
##
## 0 1
## 0 0 2
## 1 0 12par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1
## 0 1 6
## 1 0 23cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1
## 0 1 4
## 1 0 11
##
## , , = b
##
## map
## 0 1
## 0 0 2
## 1 0 12par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
multi logistic regression
# of samples is N,
# of trials of a sample i is mi,
objective variable [0,n], integer
probability of incident pi[0,1]
logit pi=log(pi/ 1-pi)=sum(bj*xji),j=[0,K],i=[1,N] (-Infinity, Infinity)
yi~B(mi, pi), # of acutual incident
odds(x)=p(x)/ 1-p(x), probablity of incident / probablity of no incident
odds ratio(x0,x1)=odds(x1)/odd(x0)
when xj -> xj +1, odds ratio -> odds ratio *exp bj ex6-3.stan
multi logistic regression
data{
int N;
int K;
array[N] int m;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
}
model{
vector[N] z=X*b;
y~binomial_logit(m,z);
}
generated quantities{
array[N] int y1;
vector[N] z1=X*b;
for(i in 1:N){
y1[i]=binomial_rng(m[i],inv_logit(z1[i]));
}
}n=30
m=floor(runif(n,1,10)) # trials are varying (1,10)
x=runif(n,-1,1)
c=sample(c('a','b'),n,replace=T)
z=x+(c=='b')
y=rbinom(n,m,1/(1+exp(-z)))
tb=tibble(x=x,c=c,y=y)
f=formula(y~x+c)
qplot(data=tb,x,y,col=c)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X,m=m)
mdl=cmdstan_model('./ex6-3.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -131.051
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 11 -87.7291 0.00014274 0.000204344 0.8569 0.8569 13
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -87.73
## b[1] 0.34
## b[2] 1.39
## b[3] 1.01
## y1[1] 5.00
## y1[2] 8.00
## y1[3] 1.00
## y1[4] 7.00
## y1[5] 0.00
## y1[6] 3.00
## y1[7] 7.00
## y1[8] 2.00
## y1[9] 1.00
## y1[10] 7.00
## y1[11] 0.00
## y1[12] 5.00
## y1[13] 2.00
## y1[14] 2.00
## y1[15] 3.00
## y1[16] 1.00
## y1[17] 3.00
## y1[18] 6.00
## y1[19] 6.00
## y1[20] 3.00
## y1[21] 4.00
## y1[22] 6.00
## y1[23] 8.00
## y1[24] 0.00
## y1[25] 1.00
## y1[26] 7.00
## y1[27] 3.00
## y1[28] 3.00
## y1[29] 6.00
## y1[30] 2.00
## z1[1] 0.47
## z1[2] 1.08
## z1[3] -0.10
## z1[4] 1.90
## z1[5] -0.09
## z1[6] 0.06
## z1[7] 2.49
## z1[8] 0.21
## z1[9] 0.36
## z1[10] 0.71
## z1[11] 0.84
## z1[12] 2.42
## z1[13] -0.70
## z1[14] 1.93
## z1[15] 0.90
## z1[16] 1.05
## z1[17] -0.20
## z1[18] 2.67
## z1[19] 0.45
## z1[20] 1.14
## z1[21] 0.08
## z1[22] 2.68
## z1[23] 1.99
## z1[24] 0.27
## z1[25] 2.16
## z1[26] 1.85
## z1[27] -0.56
## z1[28] 0.12
## z1[29] 1.02
## z1[30] 0.41mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,exc='z1',ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -89.39 -88.98 1.39 1.06 -92.11 -87.94 1.00 820 864
## b[1] 0.33 0.33 0.34 0.33 -0.24 0.88 1.00 670 595
## b[2] 1.44 1.43 0.40 0.39 0.83 2.11 1.00 1115 930
## b[3] 1.07 1.06 0.41 0.39 0.43 1.78 1.00 720 693
## y1[1] 4.26 4.00 1.38 1.48 2.00 6.00 1.00 2010 1888
## y1[2] 6.74 7.00 1.38 1.48 4.00 9.00 1.00 2086 NA
## y1[3] 3.26 3.00 1.41 1.48 1.00 6.00 1.00 1826 1859
## y1[4] 6.12 6.00 0.90 1.48 4.00 7.00 1.00 1993 NA
## y1[5] 0.47 0.00 0.50 0.00 0.00 1.00 1.00 2052 NA
## y1[6] 3.06 3.00 1.29 1.48 1.00 5.00 1.00 1915 2073
## y1[7] 6.46 7.00 0.72 0.00 5.00 7.00 1.00 1986 NA
## y1[8] 3.87 4.00 1.43 1.48 1.00 6.00 1.00 1837 1933
## y1[9] 0.59 1.00 0.49 0.00 0.00 1.00 1.00 1999 NA
## y1[10] 5.38 5.00 1.41 1.48 3.00 8.00 1.00 1965 NA
## y1[11] 0.69 1.00 0.46 0.00 0.00 1.00 1.00 1898 NA
## y1[12] 5.51 6.00 0.69 0.00 4.00 6.00 1.00 1953 NA
## y1[13] 2.62 2.00 1.46 1.48 0.00 5.00 1.00 1883 1834
## y1[14] 1.75 2.00 0.47 0.00 1.00 2.00 1.00 2025 NA
## y1[15] 4.98 5.00 1.24 1.48 3.00 7.00 1.00 2015 NA
## y1[16] 0.73 1.00 0.44 0.00 0.00 1.00 1.00 1765 NA
## y1[17] 2.63 3.00 1.28 1.48 1.00 5.00 1.00 1718 1788
## y1[18] 5.58 6.00 0.65 0.00 4.00 6.00 1.00 1836 NA
## y1[19] 5.45 5.00 1.57 1.48 3.00 8.00 1.00 1482 1745
## y1[20] 2.26 2.00 0.77 1.48 1.00 3.00 1.00 1849 NA
## y1[21] 2.60 3.00 1.18 1.48 1.00 4.00 1.00 1694 1588
## y1[22] 7.46 8.00 0.73 0.00 6.00 8.00 1.00 1661 NA
## y1[23] 7.96 8.00 1.02 1.48 6.00 9.00 1.00 1977 NA
## y1[24] 0.57 1.00 0.49 0.00 0.00 1.00 1.00 1922 NA
## y1[25] 0.89 1.00 0.31 0.00 0.00 1.00 1.00 2024 NA
## y1[26] 6.96 7.00 1.00 1.48 5.00 8.00 1.00 1842 NA
## y1[27] 1.43 1.00 1.00 1.48 0.00 3.00 1.00 1859 1781
## y1[28] 3.18 3.00 1.31 1.48 1.00 5.00 1.00 1900 1929
## y1[29] 5.93 6.00 1.28 1.48 4.00 8.00 1.00 1935 NA
## y1[30] 2.38 2.00 0.98 1.48 1.00 4.00 1.00 1965 NA
## z1[1] 0.48 0.47 0.28 0.27 0.03 0.95 1.00 1451 1682
## z1[2] 1.11 1.10 0.24 0.25 0.74 1.51 1.00 1716 1562
## z1[3] -0.14 -0.12 0.33 0.32 -0.70 0.37 1.00 782 731
## z1[4] 1.97 1.95 0.36 0.36 1.40 2.57 1.00 1444 1276
## z1[5] -0.12 -0.11 0.33 0.32 -0.69 0.38 1.00 777 724
## z1[6] 0.05 0.04 0.35 0.36 -0.52 0.64 1.00 1277 1598
## z1[7] 2.58 2.56 0.50 0.48 1.81 3.42 1.01 1303 1174
## z1[8] 0.21 0.20 0.32 0.33 -0.31 0.74 1.00 1332 1637
## z1[9] 0.34 0.35 0.34 0.33 -0.22 0.90 1.00 669 595
## z1[10] 0.73 0.73 0.25 0.25 0.33 1.16 1.00 1578 1670
## z1[11] 0.84 0.84 0.41 0.41 0.16 1.52 1.00 706 609
## z1[12] 2.51 2.49 0.48 0.47 1.77 3.32 1.01 1316 1189
## z1[13] -0.76 -0.75 0.38 0.37 -1.41 -0.16 1.00 1142 904
## z1[14] 1.99 1.97 0.36 0.36 1.42 2.60 1.00 1436 1243
## z1[15] 0.93 0.92 0.24 0.24 0.55 1.33 1.00 1672 1614
## z1[16] 1.06 1.05 0.45 0.44 0.32 1.81 1.00 756 833
## z1[17] -0.24 -0.23 0.33 0.32 -0.81 0.28 1.00 831 815
## z1[18] 2.77 2.74 0.54 0.53 1.93 3.69 1.01 1276 1158
## z1[19] 0.44 0.44 0.35 0.34 -0.15 1.00 1.00 662 600
## z1[20] 1.15 1.14 0.46 0.46 0.40 1.93 1.00 781 807
## z1[21] 0.05 0.06 0.33 0.32 -0.51 0.57 1.00 720 579
## z1[22] 2.78 2.75 0.54 0.53 1.94 3.70 1.01 1275 1158
## z1[23] 2.06 2.04 0.38 0.38 1.47 2.69 1.00 1415 1217
## z1[24] 0.27 0.26 0.31 0.31 -0.24 0.79 1.00 1355 1652
## z1[25] 2.24 2.22 0.42 0.41 1.59 2.95 1.01 1369 1165
## z1[26] 1.92 1.90 0.35 0.35 1.37 2.50 1.00 1457 1374
## z1[27] -0.61 -0.59 0.36 0.35 -1.22 -0.05 1.00 1100 863
## z1[28] 0.12 0.11 0.34 0.34 -0.43 0.68 1.00 1297 1571
## z1[29] 1.05 1.05 0.24 0.24 0.68 1.45 1.00 1702 1502
## z1[30] 0.41 0.41 0.29 0.28 -0.05 0.90 1.00 1427 1713
y1=mcmc$draws('y1')
smy=summary(y1)
table(tb$y,smy$median)##
## 0 1 2 3 4 5 6 7 8
## 0 1 0 0 0 0 0 0 0 0
## 1 0 6 0 0 0 0 0 0 0
## 2 0 0 1 3 0 0 0 0 0
## 3 0 0 2 1 1 0 0 0 0
## 4 0 0 1 1 0 1 0 0 0
## 5 0 0 0 0 1 1 1 0 0
## 6 0 0 0 0 0 1 3 1 0
## 7 0 0 0 0 0 0 0 1 2
## 8 0 0 0 0 0 0 0 1 0cat('\n')table(tb$y,smy$median,tb$c)## , , = a
##
##
## 0 1 2 3 4 5 6 7 8
## 0 1 0 0 0 0 0 0 0 0
## 1 0 4 0 0 0 0 0 0 0
## 2 0 0 0 2 0 0 0 0 0
## 3 0 0 1 1 0 0 0 0 0
## 4 0 0 1 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0 0
## 6 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0
##
## , , = b
##
##
## 0 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0 0
## 1 0 2 0 0 0 0 0 0 0
## 2 0 0 1 1 0 0 0 0 0
## 3 0 0 1 0 1 0 0 0 0
## 4 0 0 0 1 0 1 0 0 0
## 5 0 0 0 0 1 0 1 0 0
## 6 0 0 0 0 0 1 3 1 0
## 7 0 0 0 0 0 0 0 1 2
## 8 0 0 0 0 0 0 0 1 0par(mfrow=c(1,2))
hist(tb$y-smy$median,xlab='obs.-prd.',main='residual')
density(tb$y-smy$median) |> plot(xlab='obs.-prd.',main='residual')
map=c()
for(i in 1:n){
a=table(y1[,,i])
map[i]=as.integer(names(a[a==max(a)]))
}
table(tb$y,map)## map
## 0 1 2 3 4 5 6 7 8
## 0 1 0 0 0 0 0 0 0 0
## 1 0 6 0 0 0 0 0 0 0
## 2 0 0 1 3 0 0 0 0 0
## 3 0 0 1 2 1 0 0 0 0
## 4 0 0 1 1 0 1 0 0 0
## 5 0 0 0 0 1 1 1 0 0
## 6 0 0 0 0 0 0 4 1 0
## 7 0 0 0 0 0 0 0 1 2
## 8 0 0 0 0 0 0 0 1 0cat('\n')table(tb$y,map,tb$c)## , , = a
##
## map
## 0 1 2 3 4 5 6 7 8
## 0 1 0 0 0 0 0 0 0 0
## 1 0 4 0 0 0 0 0 0 0
## 2 0 0 0 2 0 0 0 0 0
## 3 0 0 0 2 0 0 0 0 0
## 4 0 0 1 0 0 0 0 0 0
## 5 0 0 0 0 0 1 0 0 0
## 6 0 0 0 0 0 0 0 0 0
## 7 0 0 0 0 0 0 0 0 0
## 8 0 0 0 0 0 0 0 0 0
##
## , , = b
##
## map
## 0 1 2 3 4 5 6 7 8
## 0 0 0 0 0 0 0 0 0 0
## 1 0 2 0 0 0 0 0 0 0
## 2 0 0 1 1 0 0 0 0 0
## 3 0 0 1 0 1 0 0 0 0
## 4 0 0 0 1 0 1 0 0 0
## 5 0 0 0 0 1 0 1 0 0
## 6 0 0 0 0 0 0 4 1 0
## 7 0 0 0 0 0 0 0 1 2
## 8 0 0 0 0 0 0 0 1 0par(mfrow=c(1,2))
hist(tb$y-map,xlab='obs.-map',main='residual')
density(tb$y-map) |> plot(xlab='obs.-map',main='residual')
gamma regression
objective variable [0,Infinity)
# of samples is N,
log (a/li)=sum(bj*xji),j=[0,K],i=[1,N]
li=a/exp(sum(bj*xji))
yi~Ga(a,li) ex6-4.stan
gamma regression
data{
int N;
int K;
vector[N] y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> a;
}
model{
vector[N] l;
for(i in 1:N){
l[i]=a/exp(X[i]*b);
}
y~gamma(a,l);
}n=20
tb=tibble(x1=runif(n,0,2),x2=runif(n,0,2),
y=rgamma(n,3,3/exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-4.stan')
mle=mdl$optimize(data=data) # it sometimes occur error and stop process## Initial log joint probability = -1039.58
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -52.3623 5.79493e-05 0.00064872 0.5011 1 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -52.36
## b[1] -1.03
## b[2] 1.59
## b[3] 1.41
## a 3.95mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.3 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -53.12 -52.74 1.60 1.26 -56.31 -51.33 1.00 618 530
## b[1] -1.00 -1.00 0.40 0.37 -1.63 -0.35 1.00 656 596
## b[2] 1.59 1.58 0.21 0.20 1.24 1.94 1.00 855 902
## b[3] 1.40 1.40 0.25 0.24 1.02 1.80 1.00 845 856
## a 3.82 3.66 1.14 1.06 2.23 5.99 1.00 945 1217
negative binomial regression
The event with probability p do not occur y times till the event occur a times
(negative binomial distribution has larger variance than poisson distribution)
y~NB(a,p), log(p)=X*b
y~NB(a,l0/(1+l0)) when y~Po(l), l~Ga(a,l0), l0=a/E[l]ex6-5.stan
negative binomial regression
data{
int N;
int K;
array[N] int y;
matrix[N,K] X;
}
parameters{
vector[K] b;
real<lower=0> a;
}
model{
a~cauchy(0,5);
y~neg_binomial_2_log(X*b,a);
}n=20
tb=tibble(x1=runif(n,-1,0),x2=runif(n,-1,0),
y=rnbinom(n,3,exp(x1+x2)))
f=formula(y~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),y=tb$y,X=X)
mdl=cmdstan_model('./ex6-5.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -169.256
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 19 -46.4117 0.000303863 6.78782e-05 1 1 23
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -46.41
## b[1] -0.57
## b[2] -2.53
## b[3] -1.22
## a 3.89mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -47.05 -46.68 1.56 1.29 -50.05 -45.27 1.00 792 1024
## b[1] -0.63 -0.62 0.51 0.49 -1.44 0.17 1.00 903 1054
## b[2] -2.60 -2.59 0.60 0.59 -3.57 -1.65 1.00 1029 1045
## b[3] -1.27 -1.27 0.55 0.52 -2.19 -0.35 1.01 1053 1023
## a 5.42 4.45 3.84 2.44 1.74 11.87 1.00 1216 859
beta regression
using prior of binomial distribution parameter p[0,1]
y~B(n,p), p~beta(a,b)
m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)
m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*sex6-6.stan
beta regression
data{
int N;
int K;
vector[N] p;
matrix[N,K] X;
}
parameters{
vector[K] be;
real<lower=0> s;
}
model{
vector[N] a;
vector[N] b;
for(i in 1:N){
a[i]=inv_logit(X[i]*be)*s;
b[i]=(1-inv_logit(X[i]*be))*s;
}
p~beta(a,b);
}n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$p)
plot(tb$x2,tb$p)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),p=tb$p,X=X)
mdl=cmdstan_model('./ex6-6.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -66.1095
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 14.9812 0.00203034 0.000370719 1 1 25
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ 14.98
## be[1] -3.28
## be[2] 6.12
## be[3] 6.63
## s 4.85mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 finished in 0.2 seconds.
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ 14.54 14.86 1.43 1.30 11.73 16.25 1.00 633 1147
## be[1] -3.35 -3.33 0.60 0.60 -4.36 -2.35 1.00 989 982
## be[2] 6.27 6.20 1.58 1.57 3.76 8.93 1.00 699 545
## be[3] 6.73 6.77 1.70 1.60 3.86 9.51 1.01 837 822
## s 4.95 4.80 1.51 1.46 2.82 7.70 1.00 1418 1481
beta binomial regression
fitting to distribution has larger variance than binomial distribution
y~betaB(n,a,b) when y~B(n,p), p~beta(a,b)
m=E[p]=a/(a+b)
s^2=V[P]=ab/(a+b)^2/(a+b+1)
m=x*be
a=ms=inv_logit(m)*s
b=(1-m)s=(1-inv_logy(m))*sex6-7.stan
beta binomial regression
data{
int N;
int K;
array[N] int y;
array[N] int n;
matrix[N,K] X;
}
parameters{
vector[K] be;
real<lower=0> s;
}
model{
vector[N] a;
vector[N] b;
for(i in 1:N){
a[i]=inv_logit(X[i]*be)*s;
b[i]=(1-inv_logit(X[i]*be))*s;
}
y~beta_binomial(n,a,b);
s~cauchy(0,5);
}n=20
tb=tibble(x1=runif(n,0,0.5),x2=runif(n,0,0.5),
p=rbeta(n,(x1+x2)*3,(1-(x1+x2))*3),
n1=floor(runif(n,5,9)),
y=rbinom(n,n1,p))
f=formula(p~x1+x2)
par(mfrow=c(1,2))
plot(tb$x1,tb$y)
plot(tb$x2,tb$y)
X=model.matrix(f,tb)
data=list(N=nrow(X),K=ncol(X),n=tb$n1,y=tb$y,X=X)
mdl=cmdstan_model('./ex6-7.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -89.2412
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 17 -60.8222 0.00598732 0.000460707 1 1 21
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -60.82
## be[1] -3.83
## be[2] 3.94
## be[3] 9.21
## s 6.08mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.2 seconds.
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 finished in 0.2 seconds.
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 finished in 0.2 seconds.
## Chain 4 finished in 0.2 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.2 seconds.
## Total execution time: 0.4 seconds.seeMCMC(mcmc,ch=F)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -60.93 -60.54 1.54 1.32 -63.98 -59.15 1.01 546 781
## be[1] -4.08 -4.08 0.78 0.76 -5.36 -2.79 1.00 856 1006
## be[2] 4.31 4.31 1.65 1.62 1.52 6.95 1.00 932 611
## be[3] 9.81 9.73 2.24 2.15 6.12 13.62 1.00 771 799
## s 49.26 10.21 455.82 7.43 3.39 84.35 1.00 915 545
